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.. . ABSTRACT r 

Smart composite structures with distribufed sensors and actuators have the capability to 
actively respond to a changing environment while offering significant weight savings and 
additional passive controllability through ply tailoring. Piezoelectric sensing and actuation 
of composite laminates is the most promising concept due to the static and dynamic control 
capabilities. Essential to the implementation of these smart composites are the development 
of accurate and efficient modeling techniques and experimental validation. This research 
addresses each of these important topics. 

A refined higher order theory is developed to model composite structures with surface 
bonded or embedded piezoelectric transducers. These transducers are used as both sensors 
and actuators for closed loop control. The theory accurately captures the transverse shear 
deformation through the thickness of the smart composite laminate while satisfying stress 
free boundary conditions on the free surfaces. The theory is extended to include the effect 
of debonding at the actuator-laminate interface. The developed analytical model is 
implemented using the finite element method utilizing an induced strain approach for 
computational efficiency. This allows general laminate geometries and boundary 
conditions to be analyzed. The state space control equations are developed to allow 
flexibility in the design of the control system. Circuit concepts are also discussed. 

Static and dynamic results of smart composite structures, obtained using the higher 
order theory, are correlated with available analytical data. Comparisons, including 
debonded laminates, are also made with a general purpose finite element code and available 
experimental data. Overall, very good agreement is observed. Convergence of the finite 
element implementation of the higher order theory is shown with exact solutions. 
Additional results demonstrate the utility of the developed theory to study piezoelectric 
actuation of composite laminates with pre-existing debonding. Significant changes in the 


ii 


modes shapes and reductions in the control authority. result due ..to partially debonded 
actuators. 

An experimental investigation addresses practical issues, such as circuit design and 
implementation, associated with piezoelectric sensing and actuation of composite laminates. 
Composite specimens with piezoelectric transducers were designed, constructed and tested 
to validate the higher order theory. These specimens were tested with various stacking 
sequences, debonding lengths and gains for both open and closed loop cases. Frequency 
changes of 15% and damping on the order of more than 20% of critical damping, via 
closed loop control, was achieved. Correlation with the higher order theory is very good. 
Debonding is shown to adversely affect the open and closed loop frequencies, damping 
ratios, settling time and control authority. 
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1. Introduction 

Structures by nature are subjected to static and dynamic forces which cause 
deformation. In the past, designers have compensated for undesirable deformation and 
vibration by adding stiffeners and passive devices which often introduce a significant 
weight penalty. The concept of smart structures has received a great deal of attention 
recently as an alternative to conventional techniques. These advanced structures can be 
designed to actively react to disturbance forces to maintain structural integrity while 
maintaining, or even improving, the level of performance. Although smart structures have 
enormous potential for an abundance of applications ranging from reducing helicopter noise 
(Anonymous; 1994), precision antenna positioning (Stevens; 1991), improving aeroelastic 
stability of fixed (Heeg; 1992) and rotory wing (Chandra and Chopra; 1993) aircraft, a 
number of basic issues still require further investigation. These issues include development 
of accurate, yet efficient, mathematical modeling techniques and experimental 
investigations. The primary focus of this dissertation is to address these key issues related 
to smart composite structures. 

The exact definition of a smart structure is rather ambiguous and deserves some 
discussion. A smart structure does not necessarily posses any inherent “intelligence”. For 
instance, the mechanical governor for a gasoline engine, which serves to automatically limit 
the speed of the engine, could be termed a smart structure. Although this device actively 
controls the operation of the engine, it is not really smart. Other useful devices have been 
invented in the past, such as the bimorph thermostat or the modern elevator, which could 
be termed smart structures. These structures actively respond to their surroundings, but do 
not posses any real intelligence. The term smart structures is a more generic term which is 
simply being used to describe the next generation of these devices aimed at enhancing 
functionality while reducing weight (Sirkis; 1996). A more precise definition has been 
suggested to characterize smart structures as load bearing structures with the ability to 



actively sense and react to its environment via onboard sensors, actuators and 
computational/control capabilities (Rogers; 1989). It must be noted that in the current 
literature, the terms smart, active, adaptive, metamorphic and intelligent are used 
interchangeably. 

Integra] components of these structures are smart materials. This class of materials can 
change their mechanical properties in response to an external stimuli. They can be used in 
sensing and/or actuation modes either independently, or as part of a larger system. Several 
different types of these materials exist, each with its own set of unique properties, and are 
described in the next section. Key issues associated with the development of smart 
structures, including mathematical modeling techniques, are discussed. 

1.1 Smart Materials 

Magnetostrictive and electrostrictive materials change shape when subjected to magnetic 
and electric fields, respectively. The most common magnetostrictive material is Terfenol-D 
which is capable of inducing strains up to 0.2%. It is named after its components terbium 
(TER), iron (FE), dysprosium (-D) and the place it was discovered, the Naval Ordinance 
Lab (NOL). Metglas is another magnetostrictive material. Although the maximum induced 
strain is much less, it can be manufactured as a foil for easy embedding in composites. An 
integrated actuation system for individual control of helicopter rotor blades using Terfenol- 
D actuators has been investigated by several researchers (Ghorayeb and Straub; 1995: 
Freidmann et ah; 1995). This class of smart materials have the advantage that wires 
connecting to the material directly are not necessary. However, the metal components used 
in these actuators result in a significant weight penalty. Other potential applications include 
torpedo control surfaces, gimballing cockpit simulators and vibration damping of optical 
benches (Card; 1992). Electrostrictive transducers have been used in a number of 
applications including adaptive optic systems, scanning tunneling microscopes and 
precision micropositioners. (Newnham and Ruschau; 1991). Both magnetostrictive and 



electrostrictive materials change shape in one direction only regardless of the polarity of the 
applied field. Therefore, a return mechanism needs to be implemented along with the 
materials themselves. This introduces additional complexities and restrictions on their 
potential for use as simple and light weight actuators. 

Another class of smart materials is electrorheological (ER) fluids which are 
characterized by a considerable variation of their rheological properties when an electric 
field is applied (Winslow; 1949). In the absence of an electric field, the strain rate of ER 
fluids is directly proportional to the applied stress, much like typical Newtonian flow 
characteristics. When exposed to an electric field, the viscosity, damping capability and 
shear strength of these materials increase significantly. These properties can also be rapidly 
altered with the application or removal of an electric field. This makes ER fluids attractive 
for providing a rapid response interface in controlled mechanical devices (Carlson et ah; 
1990). An active engine mount system has been proposed utilizing ER fluids (Card; 
1992). In this application, the viscosity of the ER fluid is adjusted to change the natural 
frequency of the system away from the dominant frequency band during service. 
Advanced helicopter rotor systems, such as soft hingeless and bearingless rotors, are 
mechanically less complex than traditional articulated rotors. However, aeroelastic 
instabilities, such as air and ground resonance, have emerged as major problems in these 
rotor systems. ER fluid based dampers are a viable candidate for active damping and 
stability augmentation (Kamath and Wereley; 1995). Cable-stayed highway bridges must 
be flexible enough to avoid damage due to seismic activity, but not so flexible as to allow 
an adverse buffet response to a wind loading. Traffic loading which generates small deck 
vibrations is also a consideration. An ER fluid damper suitable for vibration and seismic 
protection of civil structures which addresses these considerations has been designed and 
tested (Makris et ah; 1995) 



Shape memory alloys (SMA) are another type of smart material. SMAs such as Nitinol 
recover their previously undeformed shape upon heating (Funakubo; 1987). This is due to 
a change in the crystalline structure known as a reversible austenite to martensite phase 
transformation at a specific temperature. Large induced strains are possible and the 
Young’s modulus increases by almost a factor of three. SMAs have been used in 
spacecraft antenna. A wire hemisphere of the material is crumpled into a tight ball. When 
heated above 77°C, the ball opens into its original shape of a fully formed antenna 
(Schetky; 1979). In another application, the fundamental frequency of a solid rectangular 
composite beam was altered by as much as 35% using a 1.3% volume fraction of SMA 
wires (Chandra; 1993). This concept was also exploited to design a variable speed 
helicopter rotor (Chopra; 1993). Solar collectors that are focused on a central receiver are 
designed with a mechanism for defocusing the collector or deactivating it by turning it out 
of the path of the sun’s rays. This is required to avoid damaging the receiver during 
periods of inoperability. SMAs have been shown to exceed design specifications to actuate 
these solar arrays (Lobitz; 1995). Smart dental braces could also be made using SMAs 
(Newnham and Ruschau; 1991). Since SMAs are actuated by heat, the response time is 
quite slow. While electrical heating sources allow the material to respond within seconds, 
cooling often takes on the order of minutes. Although positioning and stiffening 
applications appear promising, SMAs are impractical as actuators for vibration control due 
to their slow response time. 

Currently, piezoelectric materials are the most versatile smart materials. When a 
piezoelectric material is stressed mechanically by a force, it generates and electric charge. If 
the electrodes are not short circuited, a voltage associated with the charge appears. This is 
the direct effect discovered by Jacques and Pierre Curie in 1880. A year later, Gabriel 
Lippman predicted the converse effect which was verified by the Curies. The converse 
effect occurs when a piezoelectric material is stressed electrically by a voltage which results 



in a change in the material’s dimensions (Encyclopedia Britannica; 1994). A piezoelectric 
element is therefore capable of being used as either a sensor or as an actuator when it is 
coupled to a composite substructure via surface bonding or embedding. The nature of the 
electric charge and the dimension change (positive or negative) depends on the polarity of 
the mechanical stress and the applied electric field, respectively. Consequently, 
bidirectional actuation and sensing are possible and no external return mechanism is 
necessary. This property of piezoelectric materials makes them distinct from electrostrictive 
materials which can only be actuated in one direction. Natural crystals, such as quartz, 
Rochelle salts, troumaline and lithium sulfate, were the first piezoelectric materials to be 
used. A quartz transmitter used for sonar applications appeared in 1916 (Mason; 1950). 
Early phonograph pickups were made from Rochelle salts. In the 1940s it was discovered 
that the ceramic material barium titanate could be induced to exhibit piezoelectric properties 
by exposing it to a high temperature and an electric field. Currently, most piezoelectric 
devices utilize a similar piezoceramic material, lead zirconate titanate (PZT), which was first 
used in the 1950s. The desirable properties of PZT include a high level of piezoelectric 
activity and a wide frequency range. Active flight control using PZT actuators has been 
explored by several researchers for both fixed wing aircraft (Lin and Crawley; 1995: Leeks 
and Weisshaar; 1995) and rotorcraft (Spangler and Hall; 1990: Ben-Zeev and Chopra; 
1995). Active flutter suppression using PZTs has also been investigated (Heeg; 1992). 
Vibration control using PZTs has been studied extensively (Crawley and de Luis; 1987: 
Hanagud et al; 1992: Chandrashekhara and Agarwal; 1993). PZTs have also been used in 
the emerging field of MicroElectroMechanical Systems (MEMS) for building 
microactuators (Ikuta; 1992). A video tape head positioner has been developed based on a 
bimorph PZT actuator. The nonlinear properties of some materials allows for the creation 
of tunable transducers using bias fields or forces. For instance, rubber is a highly 
nonlinear elastic medium since the material stiffens noticeably under stress. A transducer 



can be constructed which consists of alternating layers of rubber and PZT. Under low 
stress bias, the piezoelectric layers dominate at the resonant frequency of the transducer. 
Under a high stress bias, the rubber stiffens and dominates at a different resonant 
frequency. Therefore, the transducer is tunable to a specific resonant frequency using a 
stress bias (Xu; 1990). Piezoceramic materials can also be used for sensing applications. 
Accelerometers are often constructed using PZTs. A PZT sensor has been proposed to 
monitor the rate of rain fall and adjusts automobile windshield wipers to the optimal speed 
(Taguchi; 1987). PZT sensors and actuators have also been investigated for use in smart 
automobile suspensions which improve dri vability while enhancing passenger comfort 
(Tsuka and Nakomo; 1990: Thirupathi and Naganathan; 1992). 

In 1969, it was discovered that the polymer Polyvinylidene Fluoride (PVDF) can also 
develop piezoelectric properties (Kawai; 1969). The polymeric piezoelectric material PVDF 
has a variety of actuation applications including vibration control (Bailey and Hubbard; 
1985), artificial hands (Brei; 1994) and trailing edge flap actuation for rotorcraft (Seeley et 
ah; 1996). Although the compliant nature of PVDF often has certain advantages, it is not 
stiff enough to develop the force needed for most applications. Therefore, PVDF materials 
are more commonly used in sensing applications such as accelerometers (Andre et ah; 
1992) and microphones (Garner; 1977: Sullivan and Powers; 1978). PZT actuation is 
more common than PVDF due to the larger forces that can be generated resulting from the 
higher stiffness. 

1.2 Key Issues in Smart Structures 

Induced strain actuation and sensing using piezoelectric materials offers significant 
advantages, such as weight savings, over other possible actuation and sensing 
mechanisms. Therefore, piezoelectric materials are used in this dissertation. Composite 
structures, which are becoming increasingly popular due to their light weight and 
innovative design possibilities offered by ply tailoring, can be very effective as elements of 



smart structures. Many of the smart structure applications for piezoelectric materials, 
including positioning, vibration control and aeroelastic response, involve the surface 
bonding or embedding of these materials in the primary composite substructure. Several 
key issues must be addressed for the efficient implementation of these smart composite 
structures. Perhaps the most important issue is the development of accurate mathematical 
tools for the analysis of surface bonded/embedded induced strain actuators in composite 
laminates (Chopra; 1996). The mathematical analysis technique must be general in nature, 
computationally efficient and be able to include debonding between the piezoelectric 
elements and the composite substructure. Experimental investigation is also necessary to 
validate the developed model. All of these issues are addressed in this dissertation. 

1.2.1 Mathematical modeling and analysis : A key issue in the efficient implementation of 
these smart composite laminates is the development of practical mathematical modeling 
tools. There is a need for a general theory which is both accurate and computationally 
efficient and accounts for debonding at the piezoelectric actuator-composite substructure 
interface. A detailed review of the existing literature, including both analytical and finite 
element approaches, is described next. 

Bernoulli-EuJer beam analysis was used by Bailey and Hubbard (1985) to derive a 
distributed parameter control theory for active control of a cantilever beam. An uniform 
strain model was proposed by Crawley and deLuis (1987). They showed that a bonding 
layer between the piezoelectric actuator and the beam which is sufficiently thin and stiff 
becomes negligible. In a later work, Crawley and Anderson (1989) proved that the simple 
uniform strain model, which assumes a constant state of strain in the piezoelectric layers, is 
not accurate when the beam to actuator thickness ratio is less than five. They also 
discussed the nonlinear nature of piezoelectric materials at higher electric fields. The ability 
of piezoelectric actuators to produce bending and twist in a composite plate, both 
independently and simultaneously, was first investigated by Lee (1990) using classical 



laminate theory (CLT),^ His derivation of the governing equations for piezoelectric sensors 
and actuators, based on classical laminate theory, has been referred to in the literature by 
many different researchers. Wang and Rogers (1991) used Heaviside functions to model 
spatially discrete actuators embedded in composite laminates. Crawley and Lazarus (1991) 
presented an exact solution for a free isotropic plate. They also presented a Ritz solution 
for anisotropic plates with bonded actuators and more complex boundary conditions. Their 
mathematical model was based on classical laminate theory (CLT) and was verified 
experimentally. Classical approaches have also been extended to included transverse shear 
effects in anisotropic sandwich plates (Koconis et al.; 1994). All of the above analyses 
were based on the classical laminate theory which does not include transverse shear 
deformation and is therefore restricted to thin beam/plate applications. First order theory, 
also known as the Timoshenko theory for beams and the Midlin theory for plates, has been 
used to model composite structures with embedded piezoelectric actuators and sensors 
(Richard and Cudney; 1993: Tzou and Zhong; 1993). Although an improvement over 
classical theory, the first order shear deformation theory accounts for transverse shear 
deformation only in an average sense. Mitchell and Reddy (1995) recently presented a 
hybrid theory which more accurately accounts for the transverse shear strain. In their 
work, they also discuss the influence of the electric potential on a Navier type solution and 
show that it is important only in certain cases. The induced strain actuation problem for a 
beam was solved in closed form by Lin and Rogers (1993) by using approximations for the 
stress fields. The exact solution for a simply supported rectangular laminated composite 
plate, including piezoelectric actuation, has also been presented (Ray et al.; 1993). This 
work was based on an earlier elasticity solution of Pagano (1970). A power series solution 
was presented for composite cylinders with piezoelectric layers (Mitchell and Reddy; 
1995). Other analytical investigations include the optimal design of embedded piezoelectric 
actuators (Kim and Jones; 1990) and increasing the authority of actuators through discrete 



attachment techniques (Chaudhry and Rogers; 1993). The aforementioned research using 
analytical approaches provides essential physical insight into problems involving 
piezoelectric sensing and actuation. However, it is very difficult to generalize these 
analytical solutions to practical problems with realistic shapes and boundary conditions. 

The finite element method offers the flexibility to model many different types of 
structures with integrated piezoelectric materials and various boundary conditions. The 
majority of research using this approach has focused on using one and two dimensional 
beam and plate type approaches. Hanagud et al. (1992) and Hwang et al. (1993) 
developed finite element models using classical beam and plate theories, respectively. 
Viscoelastic effects have also been incorporated into a finite element model based on 
classical laminate theory (Shieh; 1993). These models are the least complex, but they do 
not account for transverse shear effects which are known to be important in the analysis of 
anisotropic composite laminates. They are also restricted to thin beams and plates. Finite 
element approaches based on first order shear deformation theory do account for transverse 
shear deformation, but only in an average sense. These approaches have been used by 
several researchers to model piezoelectric actuation of composite laminates (Shah et al; 
1993: Detwiler et al.; 1995) and have been extended to investigate vibration control 
(Chandrashekhara and Agarwal; 1993). The crude approximation of the transverse shear 
strains in the theory requires the use of shear correction factors. Furthermore, the finite 
element implementation of the first order theory is susceptible to unwarranted complexity 
and/or large inaccuracies due to spurious stiffness effects depending on the choice of 
interpolation. A very accurate approach is the layerwise theory of Reddy (Robbins and 
Reddy; 1991) which has been used to model piezoelectrically actuated beams 
(Chandrashekhara and Donthireddy; 1996) and has also been extended to included 
thermoelastic effects (Lee and Saravanos; 1995). However, the computational effort 
depends on the number of plies and can become prohibitively expensive for thick 


laminates. A global/local approach, proposed by Robbins and Reddy (1993), uses the 
concept of finite element mesh superposition in which an independent overlay mesh is 
superimposed on a global mesh to provide localized refinement in regions of interest for a 
more economical analysis. Mitchell and Reddy (1995) have developed a refined hybrid 
theory to model piezoelectric actuation and sensing in composite laminates which would be 
appropriate to implement using the finite element method. A three dimensional approach 
has also been investigated (Ha et al; 1992). An important distinction in the above analyses 
must be made regarding the electric potential. Some approaches include the electric 
potential in the total energy potential (Detwiler et al; 1995: Ha et al. 1992). The not so 
typical case of a nonuniform electric field can be modeled using this approach at the 
expense of additional degrees of freedom which are required at each of the nodes leading to 
increased CPU time. By neglecting the electric potential (Hwang et al; 1993: 
Chandrashekhara and Agarwal; 1993; Seeley and Chattopadhyay, 1996), the extra degrees 
of freedom can be neglected resulting in significant computational savings while retaining 
the ability to accurately model most practical actuator/sensor configurations. 

Relatively little attention in the literature has been paid to detailed modeling issues 
associated with adaptive composite structures, with surface bonded/embedded piezoelectric 
actuators and sensors, which include debonding. The global control capabilities of 
composite structures with piezoelectric sensing and actuation are a result of local stresses 
introduced by the piezoelectric actuators. However, the introduction of piezoelectric 
materials creates discontinuities which complicates the analysis of these smart composites. 
Debonding may also occur during the lifetime of the structure. In most of the existing 
work mentioned earlier, the actuators are assumed to be perfectly embedded or bonded to 
the primary structure . Therefore, issues associated with debonding of actuators is 
avoided. However, it has recently been shown by Seeley and Chattopadhyay (Seeley and 
Chattopadhyay; 1996) that the control authority of smart structures can be significantly 



mispredicted in the presence of debonding. The effects of these problems must be, 
investigated. Although three dimensional approaches for modeling debonding in composite 
structures (Yang and He; 1994: Whitcomb; 1989) are more accurate than two dimensional 
theories (Pavier and Clarke; 1996: Whitcomb; 1981: Kardomateas and Schmueser; 1988: 
Gummadi and Hanagud; 1995), their implementation can be very expensive for practical 
applications. The layer-wise approach (Barbero and Reddy; 1991) is an alternative since it 
is capable of modeling displacement discontinuities. However, the computational effort 
increases with the number of plies. Recently, a refined higher order theory, developed by 
Chattopadhyay and Gu (1994), was shown to be both accurate and efficient for modeling 
delamination in composite plates and shells of moderately thick construction.. This theory 
has also been shown to agree well with both elasticity solutions (Chattopadhyay and Gu; 
1996) and experimental results (Chattopadhyay and Gu; 1996). 

For the analysis of arbitrarily thick composites with surface bonded or embedded 
piezoelectric actuators and sensors, it is important to have a more effective general theory 
than is currently available. It has long been recognized that higher order laminate theories 
provide an effective solution tool for accurately and efficiently predicting the deformation 
behavior of composites laminates subjected to bending loads. However, it is difficult to 
provide a consistent displacement field which accurately satisfies the stress free boundary 
conditions at the free surfaces while maintaining continuity of strains through the thickness. 
Modeling the debonding of the piezoelectric materials has not been adequately addressed in 
the literature. To investigate this issue, a general theory is developed for the analysis of 
smart composites including the presence of debonding. The developed theory is based on 
the general higher order theory of Reddy (1990) which was implemented using the finite 
element method by Bhimaraddi and Stevens (1984) and Ren and Hinton (1986). In this 
dissertation, the higher order theory is extended to model composites with arbitrary 
thicknesses including surface bonded or embedded piezoelectric actuators and sensors. 
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The current formulation also allows for both separation and slipping which are a result of 
debonding. The theory developed here, which is implemented using the finite element 
method, is expected to more effectively model the complex stress distributions found in 
smart composite laminates with a reasonable amount of computational effort. 

1-2.2 Experimental research : Experiments are essential for validating mathematical models 
and for understanding practical aspects associated with the actual implementation of such 
structures. However, experimental studies of smart structures are reported less frequently 
in the literature than mathematical models. Bailey and Hubbard (1985) performed an 
experimental investigation to provide increased vibration damping of a flexible beam using 
a piezoelectric film (PVDF). Chopra et al. (1993: 1995: 1996) has performed several 
experiments involving piezoelectric actuation with application to rotorcraft. Crawley and 
Lazarus (1991) presented an experimental investigation of a cantilever composite plate 
utilizing piezoelectric actuation for static shape control. Torsion and bending of 
piezopolymer plates were experimentally demonstrated by Lee and Moon (1989). 
Experimental optimal control has been studied for a cantilever beam using piezoelectric 
sensors and actuators (Hanagud; 1992). Moire interferometry was used to study the 
induced strains on a plate resulting from piezoelectric actuation (Mollenhauer and Griffen; 
1994). This study concluded that numerical approaches can be used to accurately predict 
induced strains. Heeg (1992) experimentally investigated the improvement of aeroelastic 
stability for fixed wing aircraft using piezoelectric actuation. Debonding of smart 
composite laminates using piezoelectric actuation is an important issue. Yet, no 
experimental testing of composite laminates utilizing piezoelectric actuation with debonding 
has been reported in the literature. Therefore, this important topic is investigated in the 


current research. 


2. Objectives ^ 

The objective of this research is to address important issues related to smart structures. 
First, a general theory is developed to model composite structures with surface bonded or 
embedded piezoelectric transducers used as both sensors and actuators for closed loop 
control. The theory accounts for the presence of debonding. Experiments are then 
performed to validate the developed theory. 

The developed theory, which utilizes a refined higher order displacement field, 
accurately captures the transverse shear deformation through the thickness of the smart 
composite laminate while satisfying stress free boundary conditions on the free surfaces. 
The model is implemented using the finite element method to allow general laminate 
geometries and boundary conditions to be analyzed. The state space equations of motion 
are developed to allow insight into the design of the controls. The model is extended to 
incorporate the presence of debonding in the composite laminate at the interface between the 
piezoelectric actuators and the underlying substructure. The developed model is accurate, 
computationally efficient and is applicable to practical geometries. 

Extensive correlation studies are presented to demonstrate the utility of the higher order 
theory to model smart composite laminates. Comparisons, including debonded laminates, 
are made with a general purpose finite element code and available experimental and 
analytical data. An experimental investigation addresses practical issues, such as circuit 
design and implementation, associated with piezoelectric sensing and actuation of 
composite laminates. Composite specimens with piezoelectric transducers were designed, 
constructed and tested to validate the higher order theory. The composite specimens were 
tested at various stacking sequences, debonding lengths and gains for both open and closed 
loop cases. 
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Following are the principal objectives of the proposed research. 

1 . Develop a new mathematical analysis technique using the higher order theory for 
modeling smart composite laminates with embedded/surface bonded piezoelectric 
actuators and sensors. The model is accurate, general in nature and computationally 
efficient. The equations of motion include piezoelectric sensing and actuation for 
controls design. 

2. Extend the capabilities of the developed model to include debonded actuators. 

3. Correlate the developed mathematical model with published analytical and 
experimental data and results using a general puipose finite element code. 

4. Perform an experimental investigation to address practical issues such as smart 
composite laminate construction, implementation of the developed control law and 
vibration testing. Correlate the experimental results with the higher order theory. 



3. Mathematical Model Development 

A general theory is formulated to analyze composite laminates with piezoelectric 
sensing and actuation in this chapter. The theory is extended to account for debonding. 
When a piezoelectric element is stressed mechanically, it generates an electric charge. If the 
electrodes are not short circuited, a voltage associated with the charge appears. 
Conversely, when a piezoelectric element is stressed electrically by applying a voltage, its 
dimensions change. By incorporating piezoelectric materials into structures, they can be 
used as sensors and actuators to alter the static and dynamic response of the structure. The 
basis of the mathematical theory presented in this dissertation to investigate these “smart 
structures” is developed through the constitutive laws which govern piezoelectric materials 
and composite laminates. The constitutive laws are then implemented using the higher 
order laminate theory which is extended to account for debonding. The state space control 
equations are used to develop the control systems and the finite element method is used to 
implement the developed theory so that practical structures can be analyzed. 

3.1 Piezoelectric and Laminate Constitutive Relations 

The objective is to utilize the macroscopic properties of piezoelectric materials by 
integrating them as elements of smart composite structures for both sensing and actuation. 
To achieve this, it is first necessary to present the electro-mechanical constitutive 
relationships. These equations relate stress, strain, charge and electric field of a 
piezoelectric material. They are derived from the electric enthalpy density function as 
follows and were first formalized by Voigt (1928). 

H( e ij »Ej ) — ~ c jjkl £ ij£]<l ~ e ijk^i^jk — T^ijEjEj hj>k,l = 1,2,3 (3.1.1) 

where £y andE s are components of the strain tensor and electric field vector, respectively 
andCjjki, ejj k and ky are the elastic, piezoelectric and dielectric permittivity constants, 

respectively. All material constants are assumed to be isagric (measured at constant electric 



field). It must be noted that repeated indices indicate summation. The charge and stress are 
determined as follows. 
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It must be noted that there are only six unique values of the stress and strain tensors due to 
symmetry. Therefore, the following notation is used to define these quantities. 
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where x, y and z correspond to the 1, 2 and 3 directions in the tensor relationships. 
Further details of the piezoelectric constitutive relationships are found in the following 
sections and in Appendix A. 


3-1.1 Composite laminates : Consider a composite laminate which is constructed from 
layers of both orthotropic (transversely isotropic) and piezoelectric materials as shown in 
Fig 1.1. The piezoelectric coupling coefficients and the electric terms are set to zero for the 
nonpiezoelectric layers. Equation 3.1.2 is discarded and Eqn. 3.1.3 reduces to the familiar 
Hooke’s law. 

\=QijA (‘>j = 1 .2,'--,6) (3.1.5) 

where G ik are the stresses, £j k are the strains and Qy are the elastic stiffness coefficients 

defined in the material coordinate system for the k-th layer. Stresses and strains in the 
global coordinate system for the k-th layer (a ifc and £ jk ) are related as follows. 



<Ji k“^Uk 8 jk (*’J = 1,2,- * ,6) (3.1.6) 

where, in matrix form, 

Qk=T k _ 1 Q k T k 

and 

m 2 n 2 0 0 0 2mn 

n 2 m 2 0 0 0 — 2mn 

0 0 10 0 0 

0 0 0 m -n 0 

0 0 0 n m 0 

-mn -mn 0 0 0 m 2 -n 2 

where m = cos9 k , n = sin 0 k and 0 k is the angle of rotation about the Z axis between the 
global and material coordinate system for the k-th layer (Fig. 3.1). Furthermore, for 
orthotropic (transversely isotropic) plies, only five independent elastic coefficients remain 
and the transformed elastic stiffness matrix is as follows. The exact formulation of the 
elements of Qk are readily found in texts on composite materials (Vinson and Sierakowski; 
1987: Agarwal and Broutman; 1990). 

Qn Q12 0 0 0 Q !6 

Q12 Q22 00 0 q 26 

000000 
0 0 0 q 44 q 45 0 

0 0 0 Q 45 Q55 0 

Qi 6 Q26 0 0 0 q 66 





(3.1.7) 


(3.1.8) 




Fig. 3.1 Composite laminate with piezoelectric layer. 


3.1.2 Piezoelectric materials’ : Now consider the plies introduced in to the composite 
laminate which posses piezoelectric properties. Differentiating Eqns. 3.1.2 and 3.1.3 and 
utilizing the symmetry of the stress and strain tensors as before yields the following 
expressions. 

D i =e ij e j +k ik E k i = 1,2,3 (3.1.10) 

— C ij^j ~ e ki E k i — l,2,---,6 (3.1.11) 

Using the relationship 

^ij — ^ik^kj i — 1>2,3 j — 1,2,- • - ,6 (3.1.12) 


the above relationships are represented as follows. 

E i ~ ^ij c jk^k ^ k im E m * “ ^2,3 


(7j c y£j c ik^km E m 


(3.1.13) 


i = 1,2, — ,6 


(3.1.14) 
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The piezoelectric materials used in work are assumed to be either PZT (piezoceramic) 
or PVDF (piezopoly meric). Both of these materials display orthorhombic mm2 symmetry 
for which the piezoelectric coefficients contained in dy (in matrix form) are as follows. 


d - 



0 0 0 d 15 0 

0 0 d 24 0 0 

d 32 d 33 0 0 d 36 


(3.1.15) 


These piezoelectric materials are also isotropic in the context of laminate theory and are 
therefore independent of material orientation. Therefore, converting to matrix notation, 
Eqns.3.1.13 and 3.1.14 are written as follows. 
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(3.1.17) 


The above constitutive relationships are quite cumbersome in their general form. It is 
appropriate to simplify them in the context of the current research as follows. When a 
piezoelectric layer is used as a sensor, no electric field is applied (E m =0). Furthermore, the 
charge of interest (D 3 ) is normally collected on electrodes located on the upper and lower 
surfaces of the piezoelectric layer while the charge in the in-plane directions are ignored due 



to the geometry of the piezoelectric layer as shown in Fig. 1.2. This charge (D 3 ), which is 
determined for the k-th ply from Eqn. 3.1.16, is now reduced to the following form. 
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(3.1.18) 


The above equation represents the charge resulting from an applied force to the piezoelectric 
material and is referred to as the direct effect. A piezoelectric layer can also be used as an 


actuator. That is, an induced strain results from an applied electric field. Normally, the 
electric field is applied through the thickness of the piezoelectric layer used as an actuator 
and the constitutive relationship (Eqn. 3.1.18) is simplified as follows. 
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This is known as the converse effect. 
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(3.1.19) 



Fig. 3.2 Piezoelectric layer 



3.2 Analytical Development 


3.2.1 Higher order displacement field : The general higher order displacement field which 
describes deformation in the composite laminate is defined as follows. 


U(x,y,z) 

V(x,y,z) 

W(x,y,z) 


= u(x,y) + (z - c)^- — w(x,y) + cj> x (x,y) 

+ (z - c) 2 u 2 (x, y) + (z - c) 3 u 3 (x,y) 

( 3 

= v(x,y) + (z-c) w(x,y) + <j> y (x,y) 

v °y J 

+ (z - c) 2 v 2 (x, y) + (z - c) 3 v 3 (x, y) 

= w (x,y) 


(3.2.1) 


where U, V and W are the total displacements (Fig 3.3), u, v and w denote the midplane 
displacements of a point (x,y), the partial derivatives of w represent the rotations of 
normals to the midplane corresponding to the slope of the laminate and (J) x and (j) y represent 
the additional rotations due to shear deformation about the y and x axes, respectively. The 
quantities U 2 , u 3 , v 2 and v 3 represent higher order functions. The thickness coordinate, z, 
is measured from the global midplane of the laminate and c is the local midplane where c = 
0 for a laminate with no debonding present. This displacement field has the advantage of 
easily reducing to the well known classical theory if the higher order terms are eliminated. 
This particular form of the displacement field has desirable properties for the finite element 
implementation as discussed in Appendix B. In the current work, assuming that 
displacements and rotations are small, a linear relationship for the kinematic equations is 
used. 
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where E\-z 6 represent the linear strains as before. 



In the general form presented above, the higher order displacement field does not 
necessarily satisfy the requirement that the transverse shear stresses vanish on the free 
surfaces of the laminate. In addition, the displacements in the portion of the laminate where 
a debonding is assumed to be pre-existing must be independent to account for the slipping 
and separation which occurs at the debonding. These concepts are addressed in the next 
section. 


3.2.2 Refi ned displacement fields to incorporate debonding : Consider the cross section of 

a debonded smart composite shown in Fig. 3.4. The entire structure is divided into three 
regions (Fig. 3.4b) including the nondebonded region (^ u ), the region above the 

debonding (Q dl ) and the region below the debonding (^ d2 ). An interface region (S) is 

also defined where the nondebonded region and the debonded regions fo dl ,f} d2 ] 


are joined. The presence of debonding requires that the transverse shear stresses, a 4 and 
g 5 , vanish not only on the plate top and bottom surfaces, but on the debonded surfaces in 
the debonded region as well. That is, 

a 4 (x,y,a r ) = 0, a 4 (x,y,b r ) = 0 (x,y)eQ [ (r = u,dl,d2) 

a 5 (x,y.a r ) = °, 0 5 (x,y,b r ) = O (x,y)sa r (r = u,dl,d2) 


(3.2.3) 


where 


(x,y) e <T r (r - u,dl,d2) 
(x,y) e £> r (r = u,dl,d2) 
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(3.2.4 a-c) 
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For orthotropic plates, these conditions are equivalent to the requirement that the 
corresponding strains be zero on these surfaces. Therefore, 

£ 4 (x>y.a r ) = 0, s 4 (x,y,b r ) = 0 (x,y) € £2 r (r = u,dl,d2) 

e 5 (x,y,a r ) = 0, e 5 (x,y,b r ) = 0 (x,y)en r (r = u,dl,d2) 


(3.2.5) 


Three independent refined displacement fields are obtained by applying these boundary 
conditions for each iegion , r = u,dl,d2j which allows several of the higher order 

functions to be identified in terms of the lower order functions as follows. 
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where c 
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is the local midplane and d r = b r -a r is the thickness of the region 


as shown in Fig. 3.5. 
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Fig. 3.4 Smart composite cross section including debonding. 
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Fig. 3.5 Local region plate geometry. 


3.2.3 Continuity conditions : Additional boundary conditions must be imposed to ensure 
continuity of displacements at the interface of the nondebonded and debonded regions (S) 
shown in Fig. 3.6. A vector of the displacements is constructed to simplify formulation of 
the boundary conditions as follows. 

run 



(r = u,dl,d2) 


(3.2.7) 


The continuity conditions at the interface of the nondebonded and debonded regions (S), 
are imposed as follows. 

U u =U dl a dl <z<min(b u ,b dl ) (x,y)eS (3.2.8) 
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The above equations can be satisfied exactly for the classical theory since the displacement 
distribution through the thickness is linear. However, the displacement distribution using 
the refined theory is nonlinear and must be satisfied in an average sense as follows. An 
error function vector for the first of the above equations is formulated as follows. 

e = U u -U dl a dl < z < min(b dl ,b d2 J (x,y)eS (3.2.10) 


where a and b are the limits of the interval given in Eqn. 3.2.4 a-c. It is desired to minimize 
the difference between U u and U dI at each point through the thickness in S. This can be 
accomplished by first integrating the square of the error through the thickness as follows. 

E = ^J(e Te ) dz 

a (3.2.11) 


where a and b define the limits of integration through the thickness as indicated in the 
interval given in Eqn. 3.2.4 (b). These integration limits must be considered carefully 
since the presence of surface bonded actuators/sensors may change the dimensions of the 
laminate in any of the regions. It is desired to find a relationship between the independent 
displacement functions in Q dl and Q. d2 which minimizes the error in terms of the 
displacement functions in the nondebonded region to satisfy the continuity conditions. 
Therefore, derivatives of E are taken with respect to the independent functions in f2 dl and 

are set to zero as follows. 
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Taking derivatives and rearranging the above equations lead to the following relationships 
which satisfy the continuity conditions. 



3w r 

3x 

3w r 

dy 


3w u 

3x 

3w u 

3y 


+ P r 4 


r^u 
x 


+p>: 


(3.2.13 a-g) 


1>x=T>“ 

<l>y=Y r <t>; 

where r=dl and the above relationships correspond to regions H u and H dl . The constants 
a r , {3 r and Y are as follows. 


a r = (- 4a 4 - 36a 3 b - 60a 2 b 2 - 36ab 3 - 4b 4 + 36a 2 (c u ) 2 + 68ab(c u ) 2 + 36b 2 (c u ) 2 + 52a 3 c r + 
228a 2 bc r + 228ab 2 c r + 52b 3 c r - 72a 2 c u c r - 136abc ll c r - 72b 2 c u c r - 140a(c u ) 2 c r - 
140b(c u ) 2 c r - 156a 2 (c r ) 2 - 388ab(c r ) 2 - 156b 2 (c r ) 2 + 280ac u (c r ) 2 + 280bc u (c r ) 2 + 
140(c u ) 2 (c r ) 2 + 140a(c r ) 3 + 140b(c r ) 3 - 280c u (c r ) 3 - 27a 2 (d u ) 2 - 51ab(d u ) 2 - 
27b 2 (d u ) 2 + 105ac r (d u ) 2 + 105bc r (d u ) 2 - 105(c r ) 2 (d u ) 2 )/(3(9a 2 + 17ab + 9b 2 - 
35ac r - 35bc r + 35(c r ) 2 )(d u ) 2 ) 

(3 r — (18a 2 -h 34ab + 18b 2 - 35ac u - 35bc u - 35ac r - 35bc r + 70c u c r )(d r ) 2 /(2(9a 2 + 17ab + 
9b 2 - 35ac r - 35bc r + 35(c r ) 2 )(d u ) 2 ) 


Y= (- 30a 3 c u - 1 1 0a 2 bc u - 1 10ab 2 c u - 30b 3 c u + 72a 2 c u2 + 136ab(c u ) 2 + 72b 2 (c u ) 2 + 30a 3 c r 
+ 1 10a 2 bc r + 1 10ab 2 c r + 30b 3 c r + 56a 2 c u c r + 168abc u c r + 56b 2 c ll c r - 280a(c u ) 2 c r 
- 280b(c u ) 2 c r - 128a 2 (c r ) 2 - 304ab(c r ) 2 - 128b 2 (c r ) 2 + 140ac u (c r ) 2 + 140bc“(c r ) 2 + 
280(c ll ) 2 (c r ) 2 + 140a(c r ) 3 + 140b(c r ) 3 - 280c u (c r ) 3 - 18a 2 (d u ) 2 - 34ab(d u ) 2 - 
1 8b 2 (d u ) 2 + 70ac r (d u ) 2 + 70bc r (d u ) 2 - 70c r2 (d u ) 2 + 18a 2 (d') 2 + 34ab(d r ) 2 + 



1 8b 2 (d r ) 2 - 35ac ll (d r ) 2 - 35bc u (d r ) 2 - 35ac r (d r ) 2 - 35bc r (d r ) 2 + 70c u c r (d r ) 2 )/(2(9a 2 
+ 17ab + 9b 2 - 35ac r - 35bc r + 35(c') 2 )(d u ) 2 ) 

(3.2.14 a-c) 

where r=dl. Identical expressions corresponding to regions Q u and n d2 are similarly 
formulated by setting r = d2 in Eqns. 3.2. 13(a-g) It is also required that continuity of 
velocities be maintained between f2 L1 and £2 d2 . These conditions are obtained by simply 
differentiating Eqns 3.2.13 (a-g) with respect to time. Since the formulation for the 
geometric parameters is independent of time, they remain unchanged. Multiple debondings 
can be easily be incorporated into the developed theory by defining additional regions of 
debonding at arbitrary locations in the laminate. 



Fig. 3.6 Displacement distribution in cross section. 


3.3 l aminate Stress Resultants 

The laminate stress resultants, which include effects due to piezoelectric actuation, are 

formulated by integrating the stresses through the thickness as follows. 
h/2 h/2 h/2 


N: = 


Jcqdz, Mj = JzGjdz, Pj = Jz 3 Gjdz (i - 1,2,6) 


-h/2 


-h/2 


-h/2 


(3.3.1) 


h/2 

(Q].Q 2 )= j(o 5 ,0 4 )dz. 

-h/2 


h/2 


(R,.R 2 )= (a 5 ,a 4 )z 2 dz 

-h/2 


(3.3.2) 


The first three of the above terms (Nj, Mj and Pi) are the in-plane terms which are 
decomposed into two components as follows. 
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Nj = N, a - N, p 

Mj = M, a -M p (i = 1,2,6) (3.3.3 a-c) 

P = P A - P p 
1 1 L 1 A ! 

The first terms on the right hand sides of the above equations (superscript "A") represent 
the resultants from the actual strain field while the second terms (superscript "P") represent 
the resultants from the piezoelectric actuation. These in-plane terms, as well as the 
transverse shear terms (Qj, Rj i=l,2) are the laminate stress resultants and are discussed in 
the following sections. 


3.3.1 Non-piezoelectric stress resultants : It is worthwhile to write the terms from the actual 

strain field in terms of the elastic constants and strains as follows. 
h/2 

(n a ,M a ,P a ,) = f(l,z,z 3 )Q ij e j dz (i,j = 1,2,6) (3.3.4) 

-h/2 

where Qy contains the elastic constants and the strains, £j, are defined in Eqn. 3.2.2. 


Since the material properties may differ between plies, it is necessary to replace the 
continuous integrals in the above equation by the summation of integrals representing the 
contribution of each layer as follows. 

f hi- 1 


Ni 


> P i A .) = T (U^Qij^jdzl (i,j= 1,2,6) 


(3.3.5) 


hi.. 


It is convenient to write the in-plane strains in terms of their respective midplane strains and 
curvatures as follows. 

£, = £° +ZK° +Z 3 icf 
£ 2 = £? + ZK 2 +Z 3 K 2 


(3.3.6 a-c) 
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.3^2 


8* = £* + ZKf. + Z K 


where 


c o_^_ - jjx ^ 2 w 

1 3x ’ 1 3x 3x 2 


F 0_^v 
e 2 - T"» 
dy 
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5<|)y 3 2 W 

dy dy 2 


y.2 __ 4 3(|) X 


Kt — — 


3h 2 3x 
4 


o 3u 3v o 3(b v d$ v 
y = — + — , k*=— ^4- — - 
3y 3x 3y 3x 


3h 2 dy 

3 2 w 2 _ 4 

3x3y’ 6 3h 2 


3()) x , ^ 


+ 


3y 3x 


(3.3.7) 


The summation and integration in Eqn. 3.3.5 can now be carried out, allowing it to be 
written in matrix form as follows. 

N b = A b £ b (3.3.8) 


where 

N b = [n a M a P a 


(i = 1,2,6) 


(3.3.9) 


The quantity A B is the laminate stiffness matrix and £ B is the generalized strain vector. 
The laminate stiffness matrix is obtained by integrating the elastic constants through the 
thickness of the laminate, ply by ply and summing the result as follows. 


>ij] K 

D: 


sym 


E U 
F ii 
H u 


(i J = 1,2,6) 


(3.3.10) 


where 


( A ij , By , Djj , E ij , Fjj , Hjj ,) i 


hi. 


J* Qy(l,z,z 2 ,z 3 ,z 4 ,z 6 Jdz > (i,j = 1,2,6) (3.3.11 


L h k-1 
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The generalized strain vector, 
follows 


e B , is formulated using the definitions in Eqns. 3.3.7 as 
0 = 1 , 2 , 6 ) (3.3.12) 


Notice that the z-dependence is eliminated from the generalized strains by shifting this 

dependence to the formulation of the laminate stiffness matrix. 

The quantities Qj and Rj (i=l,2) arise due to transverse shear which is present in the 

laminate and are accounted for by using the higher order laminate theory. As before, these 

stress resultants are written in terms of the elastic constants and strains as follows. 
h /2 

(Q|. R i)= J(U 2 ){Q 54 £4 + Q55 £ 5} dz (3.3.13) 

-h /2 

h /2 

(Q 2 .R 2 A J* (Cz 2 ){Q 44 e 4 + Q 45 e 5 }dz (3.3.14) 

-h /2 


where contains the elastic constants and the transverse shear strains, £ 4 ^, are defined 

in Eqn. 3.2.2. Since the material properties may differ between plies, it is again necessary 
to replace the continuous integrals in the above equation by the summation of integrals 
representing the contribution of each layer as follows. 




hi. 


J (^ z2 ){Q54 e 4 + Q55 e 5 


•k-i 


(3.3.15) 


(Q 2 ,R 2 ) - = (Iz 2 ){Q44 £ 4 + Q45 £ 5} dz i 

IA-, J 


(3.3.16) 


The transverse strains are also represented by their respective midplane strains and 
curvatures. 

0 ? 2 

£4 = £4 +Z-K4 
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e 5 ~ £ 5 + z2k 5 


(3.3.17 a-b) 


and 


£4 - u 2 j , £5 - U] 


2 4 

K 4 “ ~u2 U 21 

h 


Kc = rll 


2 * U 1 1 


(3.3.18) 


This allows the transverse stress resultants to be written in matrix form as follows. 


N'p — A r p£ r p 


(3.3.19) 


where 


N T = [Qi Rif 0 = 1,2) 

The laminate stiffness matrix At is written as follows. 

>„] h f 


At = 


sym 


hi. 


(ij = 4,5) 


(3.3.20) 


(3.3.21) 


where 


(w,)=£ 


N 

k=l 


“k 

Jq,j(u 2 . 


z 4 Idz 


t n k-i 


\ (i,j = 4,5) 


(3.3.22) 


The generalized transverse strain vector is written as follows. 


N 


° K? 


0=4,5) 


(3.3.23) 


where the definitions for the midplane strains and curvatures are given in Eqns. 3.3.18. 


3.3.2 Piezoelectric stress resultants : The second terms in Eqns. 3.3.3 a-c represent the 
resultant forces due to piezoelectric actuation and are discussed next. First, these resultants 
are represented in terms of the elastic constants, piezoelectric constants and applied electric 
field as follows. 
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h/2 


(Nf.Mf.Pf^ J(l,z,z 3 )Q ij A j dz (i,j = 1,2,6) 


(3.3.24) 


-h/2 


where Qy contains the elastic constants and Aj are the induced strains due to piezoelectric 

actuation. Since the material, piezoelectric and electric properties may differ from ply to 
ply, it is necessary to integrate ply by ply and sum the contributions of each ply as before. 


W.rM-E 


N 

k=l 


hi- 


U ’ Z )Qi Jk d j k E 3 k dz 


l“ k -i 


(i J = 1,2,6) 


(3.3.25) 


where the piezoelectric constants in the k-th layer are nonzero only for the piezoelectric 
plies energized with electric field E 3) _ and expressed as follows. 


d j k ~[ d 3i d 32 0 0 0 d 36 ] k 


(3.3.26) 


It is convenient to express Eqn. 3.3.25 in terms of the applied voltage rather than the 
electric field using the following relationship. 

v k 


E 3 

J k 


(h k ~ hk-i) 


(3.3.27) 


Substituting the above expression for the electric field, Eqn. 3.3.25, can be written in 
matrix form as follows. 

hi. 





^ k = l 


'Nil 


f 

N p = 

M, p 

= 

N 

X, / 


Pf 




N 


V, 


I 


N 

k=l 


lA-i 


U' k -I 

hi. 


l n k-i 


, tQ:; d; dZ 

(h k -h k -,) Vl]k Jk 

z - — Yl q d dz 

(h k ~h k _ 3 ) VlJk Jk 

: k -Q-. d: dz 

(h k - h k _j) ^ Jk * 


(i, j = 1,2,6) 


(3.3.28) 
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It is not practical to present Eqn. 3.3.29 in a compact matrix form for the genera] case. It 
must also be noted that the terms in Eqn. 3.3.28 are nonzero only for energized 
piezoelectric layers. 

In realistic applications, piezoelectric layers are often implemented in symmetric pairs 
(Fig. 3.7). Simplification of Eqn. 3.29 is possible if the magnitude of the electric field 
applied to each piezoelectric layer is the same. Consider the common case where a 
symmetric pair of piezoelectric actuators are implemented and both the magnitude and the 
sign of the applied electric field (E3) is the same for each actuator (Fig. 3.7a). This is 
known as a unimoiph configuration and results in an in-plane force in symmetric laminates. 
Eqn. 3.3.28 is then simplified as follows. 




•< 



'[*!] 

0 0 

N p = 


Mf 


= 

0 

0 0 



K] 



0 

0 0 


A (ij = 1,2,6) 


(3.3.29) 


where 


A ; 


p y N 
ij Z^ k= i 


hi. 


J Qydz 


l n k-l 


(ij = 1,2,6) 


(3.3.30) 


The integration in the above equation is carried out only over the piezoelectric layers. The 
induced strain due to the piezoelectric actuation is as follows. 

A = [d 3l d 3 , 0 0 0 0 0 0 0] T E 3 (3.3.31) 


Bending moments are produced if one of the symmetric pair of actuators is energized 
with an equal, but opposite, electric field from the other (Fig. 3.7b). This case is termed a 
bimorph configuration and Eqn. 3.3.28 simplifies to the following. 




K] 




0 

N p = 


Mf 
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D u 



pf 




K] 


0 0 
0 

0 


A (i,j = 1,2,6) 


(3.3.32) 
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where 


h t 


D«)-X 


N 

k=l 


(z,z 3 )Qi 


dz 


(ij = 1,2,6) 


(3.3.33) 


in- 


k-l 


Again, the integration is carried out only through the piezoelectric layers and A is identical 
to Eqn. 3.3.31. It is helpful to note that the voltage, which must be applied to a 
piezoelectric layer to obtain a specified electric field, is calculated as follows 

v = E 3 t (3.3.34) 

where V is the voltage and t is the thickness of the piezoelectric layer. 




M,P 


Fig. 3. 7 Piezoelectric actuator configuration; (a) unimorph and (b) bimorph. 


3.4 Piezoelectric Sensors 

Piezoelectric layers are also used as sensors in the current research to detect the strain 
and rate of strain in the host structure. The procedure used to develop the piezoelectric 
sensor equations is based on work presented by Lee (1990) and used by Chandrashekhara 
and Agarwal (1993). It is extended here for use with the higher order laminate theory. 

3A 1-Sensor relations : When a piezoelectric layer is used as a sensor, a charge appears in 
response to the mechanical load. Only the charge in the thickness, or 3, direction (Fig. 


3.2) is of interest. The relationship for this direct effect was obtained in section 3. 1 and is 
repeated here for convenience. 


^ 3 k -dkQk e k 


where 


^k — [^31 ^32 0 0 0 d 36 ] k 


Qk 


Qn 

Ql2 

0 

0 

0 

0 


Q| 2 0 0 0 

Q22 0 0 0 

0 0 0 0 

0 0 Q 44 0 

0 0 0 q 55 

0 0 0 0 


0 

0 

0 

0 

0 


Q66 


and 


T 

£ k 



£ 2 e 3 


£ 4 e 5 



(3.4.1) 


(3.4.2) 


(3.4.3) 


(3.4.4) 


The charge enclosed by the surface S is calculated from Gauss’ Law (Lee; 1990) as 
follows 


q = 


D • ds 


s 


(3.4.5) 


where q is the charge, D is the electric displacement vector and ds is the differential area 
normal vector of S. However, this equation cannot be used directly since the results will 
be identically zero due to the fact that the charge within a dielectric is neutral. Since charge 
is built up on the surface of a piezoelectric lamina while it is under the action of an external 
force field, an equivalent circuit model is used to relate the closed-circuit charge signal 
measured from the surface electrode to the force field (Fig. 3.8). The charge built up on 
the surface due to the mechanical action is equivalent to the charge stored inside the 
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capacitor of the equivalent circuit. In order for the charge to be measured, the electric loop 
must be closed. That is, electrodes must exist on both the top and bottom surfaces of the 
piezoelectric layer so that the induced charge can move in the 3 direction. It is assumed that 
these electrodes cover the entire upper and lower surfaces, respectively, as shown in Fig. 
3.8 and are short circuited to measure the charge built up. Neglecting edge effects and 
assuming that the sensor layers are thin, the charge built up corresponding to a sensor in 
the k-th layer is approximated as follows. 


17. 



f \ 




D v. dS 

+ 

J'V lS 


► 

(3.4.6) 

vs 

j 

Z=Z k-I 

7 s J 





where S is the area of the sensor electrodes which is assumed to be the same on the top and 

bottom of the sensor layer and equal to the area of the sensor patch. Substituting the 
expression for D 3( _ into the above equation yields the following. 


1 

qk ~ 2 

r 

(l 'Qk 

• 

[ £ k-l + £ k] dS > 


l s J 


The strains are written in terms of the derivative operator matrix, H k and the vector of the 
unknown displacement functions, u as follows. 

=dlQ k |H k udA (3.4.8) 
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where 
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and 



x dx ’ y dy ’ xy dxdy 


The charge due to mechanical deformation at any given moment is related to 
which can be measured as follows. 

v k (t) = ^ 


(3.4.9) 

(3.4.10) 

(3.4.11) 

(3.4.12) 

(3.4.13) 

(3.4.14) 
a voltage 

(3.4.15) 


where Q is the capacitance of the piezoelectric layer (Fig. 3.8). This measured voltage is 
proportional to the strain in the sensor at any given time. The current produced is obtained 
by differentiating the charge with respect to time as follows. 


i k( t) = fM 

kW dt 
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(3.4.16) 


The current is proportional to the strain rate. 


3 


1 



piezoelectric sensor equivalent circuit 

Fig. 3.8 Equivalent circuit for piezoelectric sensor in k-th layer. 


Sensor circuits : The ultimate goal of implementing piezoelectric sensors is to obtain a 
signal, in the form of a voltage, which is proportional to either the strain or strain rate 
present in the structure. This is accomplished by connecting the piezoelectric sensors to 
op-amp circuits which can then be coupled to either a data acquisition system or a feedback 
control loop. Strain measurements are accomplished through the use of a charge amplifier 
shown in Fig. 3.9. The output voltage, vq is calculated as follows. 

_ q(t) 


v o(t) 


c, 


(3.4.17) 


The circuit must be designed around the lower and upper cutoff frequencies given as 
follows. 


f ] = 1 

1 cp 


27tR f Cj 


(3.4.18) 


f 2 


cp 2 tiR,C 


(3.4.19) 
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1 2 

where f cp and f cp are the lower and upper cutoff frequencies, respectively. The feedback 

capacitor, Cf, is chosen to match the capacitance of the piezoelectric sensor, C. Once the 
cutoff frequencies are determined, R\ and Rf are obtained from Eqns 3.4.18 and 3.4.19. 


Cf 



Fig. 3.9 Piezoelectric sensor circuit for strain measurement. 


A transimpedance amplifier (Fig. 3.10) is used to convert the sensor current, which is 
proportional to the strain rate, to a voltage which can then be used by the measurement 
instrumentation. The output voltage, v 0 , is proportional to the sensor current, i, as 
follows. 

v o(0 = -R 2 i(t) (3.4.20) 

v 0 = -R 2 i 


(3.4.21) 



equivalent piezo sensor circuit 


r 2 



Fig. 3. 10 Piezoelectric sensor circuit for strain rate measurement. 


3.5 Equilibrium Equations and Boundary Conditions 

A variational approach using Hamilton’s principle is used to derive the equilibrium 

equations and boundary conditions for a nondebonded plate including piezoelectric 

actuation. The formulation can easily be extended to debonded plates by specializing the 
following analysis to each region (n r ,r = u,dl,d2). In the absence of any 

nonconservative forces, Hamilton’s principle is stated as follows. 

8j[T-(V + U)]dt = 0 (3.5.1) 

h 

where T is the kinetic energy, V is the potential energy due to external forces, U is the 
strain energy and t] and t 2 are the initial and final times, respectively The detailed 
derivation of these quantities are presented in the following sections. 

3.5.1 Potential energy formulation : The elastic strain energy, U, is calculated by 
integrating the strain energy density, Uq, over the volume of the laminate as follows. 
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h/2 


u = | j Ju 0 dAdz 


-h/2 A 
h/2 


\ j J[ a l £| +° 2 e 2 +° 3 e 3 + ° 4 e 4 +a 5 e 5 +<7 6 £ 6 ]dAdz 


-h/2 A 


The variation of the potential energy is obtained as follows 
h/2 

5U = J J(a 1 5e 1 +a 2 8£ 2 +a 4 5£ 4 +a 5 6e 5 +a 6 5£ 6 )dAdz 


-h/2 A 


(3.5.2) 


(3.5.3) 


It is noted that the transverse normal strain associated with the assumed displacement field 
(Eqn. 3.5.2) is zero. Hence, the admissible virtual strain is also zero, making the term 
g 3 5£ 3 = 0 in Eqn. 3.5.3. On the other hand, o 3 is assumed to be negligible by using the 

plane stress assumption. Integration through the thickness and substitution of the 
expressions for the variations in the strains, obtained from Eqn. 3.2.2, yields the following 
expression. 
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(3.5.4) 



where N,, Mj and Pj (i=l,2,6) and Qj and R\ (i=l,2) are the laminate stress resultants 
which are defined in Eqns. 3.3.3 

The potential energy due to an applied distributed load is formulated as follows. 

V = Jp(x,y)U 3 dA (3.5.5) 

A 

where p(x,y) is the distributed load. Upon substituting the appropriate expression for U 3 , 
taking the variation of V yields the following expression. 


5V = 


p(x,y) 6 wdA 


(3.5.6) 


3.5.2 Kinetic energy formulation : The kinetic energy, T, is calculated by first considering 
the kinetic energy of a system of particles relative to an inertial frame of reference (Fig. 
3.11) 

T = “X mi ' 7i#?i (3-5.7) 

i=l 


where v ( is the absolute velocity of the i-th particle. 
Vi = Vj +(£2; Xl-j) 


(3.5.8) 
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Fig. 3.11 Particle relative to inertial frame of reference. 

In the context of laminate theory the rotational velocity term x r ; ) is much smaller than 
the transitional velocity term (vj). Therefore, it is appropriate to retain only Vj. 

1 

T = -2, ra i v i* v i (3-5.9) 

i=l 

Summation of all of the infinitesimally small particles throughout the laminate allows the 
sum to be represented as an integral over the volume of the laminate as follows. 

T= IJ (v * v)dV (3 ' 5 - 10) 

v 

where the transitional velocity, v, is calculated by taking the time derivative of the refined 

displacement field. Carrying through the dot product in Eqn. 3.5.10 , substituting in terms 

from the displacement field and taking the variation of T yields the following expression. 
h/2 

-4 J i I 

-h/2 A 



( 3 . 5 . 13 ) 



3.5.3 Hamilton’s principle : The variation of the integral in Eqn 3.5.1 is performed as 

follows. 

h 

f[ST - (5V + SU)]dt = 0 (3.5.14) 

ti 

The formulations for the variation of the kinetic energy (ST), the potential energy due to an 
applied load (5V) and the variation of the strain energy (8U) are now substituted into Eqn. 
3.5.14. Further integration by parts using Green’s theorem yields the following set of 
governing equations and boundary conditions. 
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The boundary conditions are as follows, 
specify: 
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where n x = cos(c|)), n y = sin((j)) and <(> represents the angle tangent to the boundary. 
Although the relevant equilibrium equations and boundary conditions are quite lengthy, 
important physical interpretation of the contribution of the higher order terms can be 
gained. As a check, when the higher order terms are set to zero, the above equations 
reduce exactly to those obtained using the classical theory. It must be noted that in the 
derivation of the equations of motion, the influence of the electric potential is assumed to be 
small and is therefore neglected. In the literature, the treatment of the piezoelectric effects 
in this manner is referred to as an induced strain formulation (Mitchell and Reddy; 1995). 

3.6 Finite Element Implementation 

Application of the higher order theory, including piezoelectric sensing and actuation, to 
realistic problems requires that complex geometries and boundaries be incorporated into the 
analysis. The logical choice for the solution technique is the finite element method which 



easily accommodates material discontinuities and is well suited for numerical 
implementation and solution. The following sections contain details of the finite element 
formulation. 


3.6.1 Equations of motion : The displacements in matrix form are presented as follows. 



(3.6.1) 


where the subscript “r” corresponds to the nondebonded and debonded regions 
(r = u,dl,d2,...). In the following development, it is assumed that specialization to each 
region is performed when necessary. Therefore, the subscript “r” is dropped. The global 
displacements are represented in matrix form as follows. 

U = GL g u (3.6.2) 


where u is a vector of the five unknown displacement functions, 
r -tT 


u — 


u 


<t>x 


V (j)y W 


(3.6.3) 


The quantity L G is a derivative operator matrix and G contains the z dependence of the 
displacement field. These matrices are defined as follows. 


L g = 


1 0 0 0 0 

0 10 0 0 
0 0 10 0 
0 0 0 10 
0 0 0 0 1 

0 0 0 0 d x 

0 0 0 0 d y 


(3.6.4) 
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1 



z — z 3 — — y 0 0 0 -z 0 

3d 2 

0 1 z - z 3 —if 0 0 -z 

3d 2 

0 0 0 10 0 


and 



^ , _ 9 2 , _ , _ a 2 

dy ’ xx Ox 2 ’ yy dy 2 ’ xy dxdy 


(3.6.5) 


(3.6.6) 


z = z-c (3.6.7) 

where c is the local midplane and d is the local thickness. Theses quantities vary from 
region to region. The generalized in-plane and transverse strains (Eqns. 3.3.6 and 3.3.17) 
are written in matrix form as follows 


£ b _ L b u (3.6.8) 

£t = L t u (3.6.9) 

where L B and L T are derivative operator matrices. The derivative operator matrices are 
defined as follows. 
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(3.6.10) 


L*t — 


0 0 
1 0 
0 0 


0 — T o 


1 01 
0 0 

4 0 

d 2 

0 0 


The kinetic energy is presented in matrix form using Eqn. 3.5. 10 as follows. 

T = — I" pU T UdV 
2 * 

V 


iJpii T L T G GL G udA 


where 


j 


G = I G T Gdz 


Using Eqns. 3.6.8 and 3.6.9, the strain energy is written as follows 


1 


J £ B^BdA + J l 


U = l\) e B N B dA + J 4 N T dA 

I A A 


J* EgAgSgdA — J* £ B N p dA + J* £^Ay£ydA 


IA 


Also, the potential energy due to an applied transverse load, p(x,y), is as follows. 


V = J U 1 FdA 

A 


rTi 


where 


49 

(3.6.1 1) 

(3.6.12) 

(3.6.13) 

(3.6.14) 

(3.6.15) 

(3.6.16) 

(3.6.17) 


F = [0 0 p(x,y)] 


(3.6.18) 
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Each of the five unknown functions are represented by their corresponding elemental 


functions u e (x,y) which are interpolated as follows. 
N„ 

u e (x,y) = ^^Nf(x, y )wf 


(3.6.19) 


i = l 


where N n is the number of nodes, N? are the interpolation functions and the superscript e 
denotes the corresponding parameter at the element level. The quantities wf are the nodal 
degrees of freedom defined as follows. 


wf — 


u i <i>xi v i $ 


w ■ 
yi w i 


3wj 3wj 


~|T 


3x Oy 


(3.6.20) 


Bilinear shape functions are used for the first four unknowns while 12 term cubic 
polynomials are used for the transverse displacements (w). The derivative terms are 
included as degrees of freedom since C 1 continuity is required for the transverse 
displacement field . The resulting four noded rectangular elements are nonconforming for 
computational efficiency and contain 28 degrees of freedom each. 

The elemental (U e ) displacements are now represented as follows 


U e = GBq^ 


, e 


(3.6.21) 


where 


Br. = L n N e 


(3.6.22) 


The elemental generalized in-plane and transverse strains are represented as follows 


8 e B =B C B W e 


(3.6.23) 


Ej = B^w e 


(3.6.24) 


where 


Bb - Lb 


N e 


(3.6.25) 
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B t = l t n ° ~ ( 3>6 26 ) 


The elemental expressions for the displacements and strains are now substituted into the 
formulation for the kinetic, strain and potential energies. (Eqns. 3.6.12,15,17) Taking the 
variation of each quantity yields the following expressions. 


8T = Sw eT M c w e 

(3.6.27) 

SK = 8w cT K e w c - 6 F; t 

(3.6.28) 

8V = 8w cT F c 

(3.6.29) 

where 


M e = jBg'GB^dA' 

(3.6.30) 

Ae 


K e = I B| T A E B|dA e + f B^- t A T BjdA c 

J J 

(3.6.31) 

A e A e 


F e = f N eT p e (x,y)dA e 

(3.6.32) 

Ae 


Ff, = J N eT N p dA e 

(3.6.33) 


A e 


The quantities M e , K e and F e represent the elemental mass and stiffness matrix and force 
vector due to a distributed load, respectively. The quantity Fp is the force vector due to the 
piezoelectric actuation. The laminate stiffness matrices (Ap, Ay) are integrated analytically 
through the thickness of the laminate, ply by ply and the finite element matrices are 
assembled using exact or full numerical integration. It must be noted that reduced 
integration and shear correction factors are not needed in the current analysis. According to 



the discretized form of Hamilton’s principle where the above quantities are summed over all 
of the elements, the following must be true 

8n = J' 2 ^[ST e - 8U e + sv e 




t = 0 


(3.6.34) 


c- 


where t\ and tj are the initial and final times, respectively and N e is the number of 
elements. Integration by parts and consideration of the arbitrary nature of the variation 5w e 
leads to a linear set of equations which are solved for the nodal displacements w. 

Mw + Kw = F + F P (3.6.35) 

where the quantities M, K , F and w denote the global mass and stiffness matrices, the 
force vector due to a distributed load and the nodal displacement vector, respectively. The 
quantity Fp is the force vector due to the piezoelectric actuation. 


3.6.2 Sensor equations : The charge from the piezoelectric sensors is determined 

independent of the finite element equations in post computation calculations since the 
electric displacements are not considered as additional degrees of freedom. This leads to a 
significant amount of computational savings and is consistent with the current approach 
which is an induced strain formulation. The short circuit voltage in the k-th layer is related 
to the charge and the displacements as follows. 

V K =^r = ^- d kQkJ H k udA (3.6.36) 

k k 


where V k is the voltage, q k is the charge and C k is the capacitance of the k-th layer. The 
displacement vector u is represented in terms of the shape functions N e and nodal 
displacements, w e , as before. Summing over all elements and plies leads to the following 


expression. 
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N, N 


v = IIf^ HCw = 

'-k 


(3.6.37) 


e=l k = l 


where 

H e 


■; 


HN e dA c 


(3.6.38) 


It is important to note that only the piezoelectric plies which are used in a sensor mode are 

to be included in the above sum. The current is determined as follows, 
dq 

(3.6.39) 


i(t) = 


dt 


Summing over all plies and differentiating yields the following expression. 


(3.6.40) 


e=l k=l 


3.7 Implementation of Continuity Conditions 

One the finite element model has been constructed, it is necessary to implement the 

continuity conditions to ensure continuity of displacements at the interface of the 

nondebonded and debonded regions (S). For simplicity, consider the simple case where 
the displacements in a nondebonded region, (^ u ), must be identical to the displacements in 

the debonded region, at the interface between these two regions, (S) (Fig. 3.12). 

Since the displacements are represented in terms of the nodal quantities in the finite element 
implementation, this is accomplished by applying the continuity conditions developed in 
Section 3.2.3 on the unknown displacement quantities associated with the nodes contained 
in S. These constraints could be applied directly using the Lagrange multiplier technique. 
However, this leads to a nonsymmetric set of equations which has very undesirable 
consequences in the solution sequence of the finite element model. A more efficient 
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approach which retains symmetry of the equations of motion needs to be developed. This 
is outlined below. 



Fig. 3.12 Finite element discretization of debonding. 


The continuity conditions presented in Eqns. 3.2.13 a-g, between the nondebonded 
(n u j and debonded (<T d ) regions, are applied to the finite element degrees of freedom at 

the interface of the nondebonded and debonded regions (S) by first presenting these 
discretized conditions in matrix form as follows. 
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(3.7.2) 


W = I u u u d 


(3.7.3) 


u = 


u u <$>" v u <J)“ w l 


3w u 3w u 
3x dy 


(3.7.4) 
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(3.7.5) 



c = c 


^ _ 55 
(3.7.6) 


Quantities with superscripts u and d refer to the corresponding quantities contained in S as 
well as the nondebonded (e> u ) and debonded (^ d ), regions, respectively. The 

formulation for a r (3 r and y are presented in Eqns. 3.2.14 a-c. It must be noted that to 
obtain Eqn. 3.7.1, Eqns. 3.2.13 (d-g) have been multiplied by the thickness of the 
debonded layer, d r , for consistency of units. Continuity of velocities must also be 
maintained. Therefore, Eqn. 3.7.2 is differentiated with respect to time to yield the 
following expression 

Rw = 0 (3.7.7) 


Next, the discretized potential energy is reformulated as follows. 


y - i-^- w ] Kw + w t F + w T F p )-— pjW T R T Rw 


(3.7.8) 


where the first term in parentheses on the right hand side of Eqn. 3.7.8 corresponds to the 
actual potential energy while the last term corresponds to the penalty term related to the 
continuity constraints which has been introduced. Minimization of this penalty term leads 
to satisfaction of the continuity conditions. Similarly, the kinetic energy is reformulated as 
follows. 


T = 


— w t Mw + — p 7 w T R 1 Rw 
2 2 2 


(3.7.9) 


where the first term on the right hand side represents the actual kinetic energy while the 
second term is the penalty term related to the continuity constraints. Again, minimization of 
this penalty term requires that the continuity conditions be satisfied. The scalar penalty 
factors, pi and p 2 , are chosen to be on the order of the 1-norm of the stiffness and mass 
matrices, respectively. According to Lagrange’s method for a discrete system, 


d far ' 


3T 9L nc 

^r Qk 


(3.7.10) 
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where are the generalized coordinates and QjJ c are the generalized forces. The 
following symmetric set of augmented equations of motion is obtained ^with Q£ c = 0^ 

[M + p 2 P]w + [K + pjP]w = F + F P (3.7.11) 

where P is the penalty matrix corresponding to the continuity constraints be satisfied. It is 
obtained first in terms of the nodes contained in (S) as follows. 

P = R T R (3.7.12) 

The penalty matrix P is then expanded to correspond to the global degrees of freedom (u) 
and has the following form. 
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(3.7.15) 


(3.7.16) 


Rji = Rj 


(3.7.17) 


3.8 State Space Controls Analysis: 

The discretized linear equations of motion to model piezoelectric sensing and actuation 
of composite laminates, including debonded, are stated as follows. 

M * w + Cw + K * w = F + F P (3.8.1) 

where 

M* = [M + p|P] (3.8.2) 


K* = [K + p 2 P] (3.8.3) 

The above matrices are determined from the finite element implementation of the higher 
order theory. It must be noted that the viscous damping matrix C has also been introduced. 
This damping term arises due to internal friction. A simple damping model is used for 
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computational efficiency. The static displacements, w, are obtained by solving the 
following set of linear equations. 

K * w - F + F P (3.8.4) 

The undamped open loop frequencies and mode shapes are calculated by solving the 
following eigenvalue problem. 

(-cofM*+K*)(|) i =0 (3.8.5) 

where CDj are the undamped natural frequencies and (j)j are the undamped mode shapes for 
the i-th mode. The eigenvectors (mode shapes) are combined into a single modal matrix as 
follows. 

<I , = [4>1 <t>2 ••• <Cn] (3.8.6) 

where N < NDOF, the total number of degrees of freedom contained in Eqn. 3.8.1. The 
mode shapes are orthogonal and are mass orthonormalized as follows. 

=1 (3.8.7) 

Classical damping, also known modal damping, is assumed (Meirovitch; 1990). Using 
this approach, the undamped and damped mode shapes are identical and remain orthogonal. 
The damping matrix also has the following property 

O t CO = Z (3.8.8) 


where 




Z = 


2 ^ 2©2 


(3.8.9) 


2?n©n_ 


and are the damping ratios for the i-th mode. The original discretized degrees of 
freedom, w, are represented in term of the mode shapes and modal participation factors as 
follows 


w = ®q 


(3.8.10) 


where q is a vector of the modal participation factors. 

"qi ’ 

Q2 


(3.8.11) 


On 


By pre multiplying the original equations of motion by O t , the modal equations of motion 
are written as follows. 

q + Zq + Aq = O t [f + F p ] (3.8.12) 


where 


co 


2 
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A = <J> T K<t> = 




(3.8.13) 


In practice, potentially thousands of degrees of freedom contained in Eqn. 3.8.1 are 
accurately represented using just a few modes in Eqn. 3.8.12. The response of the 
piezoelectric sensors is now represented as follows. 

v = H v Oq (3.8.14) 

i = H c Oq (3.8.15) 


where v is the vector of one or more sensor voltages and H v is the global sensor matrix 
which relates the sensor voltages to the displacements. Similarly, i contains the sensor 


currents which are obtained using the global matrix H c which relates the sensor currents to 
the velocities. 

In state space, the above equations are represented in standard form as follows. 


d_ q 

dt q 



+ Bu 


(3.8.16) 


y = c 


q 

q 


(3.8.17) 


where A is the plant matrix, B is the control matrix and C is the observer matrix. The 
quantity q contains the modal participation factors, u contains the control inputs and y 
contains the sensor outputs. 

These quantities are expressed as follows. 
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(3.8.18) 
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B <]) t F ® t f' p 

(3.8.19) 

C = [H V 4> H C <D] 

(3.8.20) 
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(3.8.21) 


1 



(3.8.22) 


The first element of the control vector, u, is simply the number one which corresponds to 

the external disturbance force vector, F, contained in the control matrix (B). It must be 
noted that the quantity F p contained in B are the control forces per unit volt due to each 

piezoelectric actuator. The input voltages to the actuators are contained in u ( vj , v^, etc.) 



The voltage inputs are related to the sensor voltage and current outputs for feedback control 
through a gain matrix G. 

u = Gy (3.8.23) 

The standard state space notation used here to represent the equations of motion and the 
feedback control allows the effects of piezoelectric sensing and actuation to be used in the 
design of a multi input multi output (MIMO) control system 

3.9 Implementation of the Finite Element Model: 

3.9.1 Laminate discretization : A few words regarding the discretization of the 
composite laminate incorporating piezoelectric sensors and actuators for use with the finite 
element model are in order. First, boundaries of the piezoelectric sensors and actuators 
must correspond to element boundaries as shown in Fig. 3.13. Although the piezoelectric 
sensors and actuators may occupy one or more elements, a single element either does or 
does not contain one ore more piezoelectric layers. 

Although small meshes can be formulated by hand, the generation of large meshes with 
multiple sensors and actuators quickly becomes a daunting task. An automatic mesh 
generation procedure to incorporate piezoelectric transducers is essential, although it has 
never been discussed in the literature. Therefore, a new mesh generation procedure is 
developed which is described in Appendix C. 
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(a) valid discretization (b) invalid discretization 

Fig 3. 13 Finite element discretization of piezoelectric sensor/actuator. 


The next issue involves the modeling of debonding. The higher order theory and 
corresponding finite element implementation are essentially two dimensional in nature. 
Nodes are defined at the laminate midplane. All of the unknown quantities are calculated in 
terms of these in-plane nodal values from the finite element solution. However, a 
debonding is assumed to occur between two plies at an arbitrary location through the 
thickness as well as the in-plane directions. Therefore, it is necessary to define the nodal 
points in terms of the thickness (z) direction as well. The regions both above and below 
the debonding are also discretized at their corresponding midplanes as shown in Fig. 3.14. 
The continuity conditions, as described above, are imposed on the nodes at the interface 
between the nondebonded and debonded region (S). 





Fig 3.14 Finite element discretization of debonding 

3.9.2 Computational considerations : The mass and stiffness matrices which are used in the 
finite element solution (FEM) of the static and dynamic equations of motion (M* and K*) 
are potentially quite large. Therefore, solution of the equations of motion may become 
computationally prohibitive. Some consideration of the form of these matrices leads to a 
considerable savings in computational effort. These matrices are in general both symmetric 
and sparse. Storage requirements of the entire matrices (dense storage) increases as N 2 , 
where N is the size of M* and K*. Storage of only the symmetric terms gives some 
degree of savings. However, storage requirements of only the nonzero symmetric terms 
(sparse storage) increases at a much more manageable linear rate as shown in Fig. 3.15 and 
gives a significant savings over other types of storage schemes. 

The displacement vector is the unknown quantity in the static equations of motion (Eqn. 
3.8.4). The solution to these linear equations is obtained using any one of several standard 
iterative technique for solving large sparse symmetric linear systems. The Jacobi 
Conjugate Gradient method is selected for this research (Kincaid et ah; 1996). The 
generalized eigenvalue problem found in Eqn. 3.8.5 is solved using the Lanczos method 
which is most effective for finding a few eigenvalues and eigenvectors of the large sparse 
symmetric generalized eigenproblem (Jones and Patrick; 1996). 





- = •- 4. Verification Studies 

The higher order theory, as implemented using the finite element method in this 
research, must be correlated with other approaches to ensure its accuracy. Therefore, the 
purpose of this chapter is to provide validation studies to access the validity of the higher 
order theory for modeling piezoelectric actuation of composite laminates with debonding. 
Verification is established for the fundamental cases of plates with no debonding or 
piezoelectric actuation before additional complexity is added. Then, major features of the 
new theory, such as piezoelectric actuation and debonding, are introduced. Exact 
approaches provide excellent benchmarks to access the validity of the higher order theory 
when available. Experimental data is also very useful for comparison. Other approaches, 
such as the Ritz method and commercial finite element codes, such as NASTRAN, provide 
useful verification tools as well. Standard isotropic plate problems are used in the first part 
of the validation. Material orthotropy is then introduced. Further comparisons are then 
presented which include piezoelectric actuation and finally debonding. Comparison of the 
current higher order theory with other models is also discussed. 

4.1 Fundamental Correlation of Isotropic Plates 

4.1.1 Simply supported isotropic plates : In this section, the developed higher order theory 
(HOT), implemented with the finite element method, is compared with exact solutions for 
the classical laminate theory (CLT), first order shear deformation theory (FSDT) and exact 
elasticity solutions. The exact solutions for the CLT and FSDT are obtained according to 
standard techniques (Reddy, 1984). A 4x4 quarter plate finite element mesh is used to 
generate the static HOT solution while and 8x8 mesh for the full plate is used to obtain the 
plate natural frequencies so that anti symmetric modes will not be missed. The test article is 
an isotropic, simply supported plate with material properties E = 74.5 GPa and v = 0.3 and 
an applied uniform distributed load shown in Fig 4.1. The plate is analyzed for several 
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different widths and thicknesses using CLT, FSDT and IIOT.'"' The static results are 

recorded in Table 4.1 where the nondimensional displacement (w*) is presented at the 
middle of the plate and is formulated as follows. 

„ w(a/2,b/2)h 3 E 

w* = 4 

q ° a (4.1.1) 

where q=qo is the uniform distributed load. Natural frequencies are presented in Table 4.2 
where the first 10 nonrepeated frequencies are presented. The frequencies are normalized 
as follows. 

a 2 fcT 

03 i * = C 0 i — i = 1,2,3,... (4.1.2) 

h v G 

All theories agree quite well for the static results for thin plates as expected (a/h=100). 
However, significant deviations between CLT and the other theories occur for thicker 
plates at lower a/h ratios, even for this simple isotropic example. This is because classical 
theory does not take into account transverse shear deformation and is therefore not valid for 
thick plates. The FSDT is able to account for a constant transverse shear stress through the 
thickness. The HOT, which has the same number of unknowns as the FSDT, allows for a 
quadratic variation in transverse shear deformation and predicts a significantly larger 
displacement compared to CLT, but agrees with FSDT. The HOT consistently gives a 
slightly higher deflection, but never deviates more than 2% from the FSDT. The HOT is 
also most accurate for a/b=l since the elements are square. At other a/b ratios, the aspect 
ratio of each element is no longer equal to one and the solution is slightly less accurate. 
These results indicate that the HOT agrees with exact solutions for CLT and FSDT for thin 
plates as it should, but deviates from the CLT while still agreeing with the FSDT for thick 
plates. 
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Fig. 4.1 Simply supported isotropic plate with distributed load. 


Table 4.1 Comparison of normalized center deflections (w*) of simply supported isotropic 
plates under uniform distributed load. 


b/a 

a/li 

CLT 

w* 

FSDT 

HOT 


100 

0.0444 

0.0444 

0.0448 

1 

10 

0.0444 

0.0467 

0.0471 


5 

0.0444 

0.0536 

0.0541 


100 

0.1106 

0.1106 

0.1119 

2 

10 

0.1106 

0.1142 

0.1155 


5 

0.1106 

0.1248 

0.1262 


The normalized natural frequencies are presented in Table 4.2. for the CLT, FSDT, the 
exact elasticity solution and the current HOT using an 8x8 mesh. The exact elasticity 
solution is obtained by formulating the equilibrium equations, where no assumptions are 
made for the displacements, and solving using Fourier series (Srinivas et al.: 1970). The 
CLT consistently over predicts the frequencies for this moderately thick plate (L/h=10) 
since it neglects transverse shear deformation. The HOT agrees well with both the FSDT 
and the exact elasticity solution. The deviation increases with mode number since the mesh 
for representing the higher modes is somewhat crude. Errors are within acceptable limits 
and never exceed 5% for the modes presented. 








_ 68 
Table 4- 2 Comparison of normalized natural frequencies of simply supported isotropic 
plates (a/h=10, b/a=l) 


m 

n 

CLT 

FSDT 

Exact 3 

HOT (8x8) 

1 

1 

0.0963 

0.0930 

0.0932 

0.0922 

2 

1 

0.2408 

0.2219 

0.2260 

0.2192 

2 

2 

0.3853 

0.3406 

0.3421 

0.3304 

1 

3 

0.4816 

0.4149 

0.4171 

0.4100 

2 

3 

0.6261 

0.5206 

0.5239 

0.5004 

1 

4 

0.8187 

0.6520 

- 

0.6469 

3 

3 

0.8669 

0.6834 

0.6889 

0.7169 

2 

4 

0.9632 

0.7446 

0.7511 

0.7375 

3 

4 

1 .2040 

0.8896 

- 

0.8611 

1 

5 

1.2521 

0.9174 

0.9268 

0.9207 


a) From Srinivas et al. (1970) 


4.1.2 Static and dynamic convergence evaluation : In this section, convergence of the 

finite element model developed for the higher order theory (HOT) is investigated for both 
static and dynamic analyses. The test article is a simply supported isotropic plate with a 
constant distributed load identical to the one used in the previous section. For both static 
and dynamic cases, the HOT is compared to the corresponding exact solution for first order 
shear deformation theory (FSDT). 

Convergence for the static case is investigated by examining a quarter plate model 
where the aspect ratio b/a corresponds to the full plate. In Fig. 4.2, the HOT normalized 
displacement results are presented for six meshes ranging from lxl to 6x6 for several 
different thicknesses with b/a=l. The displacement is normalized in the same manner as 
the previous section. These results are compared to the corresponding exact solutions 
using FSDT. In all cases, convergence is very rapid. The 4x4 mesh yields acceptable 
results with very little CPU time required. The 6x6 mesh is slightly more expensive, but 
yields displacements which are within 1% of the exact solution. The HOT also converges 






0.17 — i 


— d — HOT (a/h -TOO) 



1 i I 1 1 1 1 

0 lxl 2x2 3x3 4x4 5x5 &6 


Mesh size (quarter plate, b/a = 2) 

Fig. 4.3 Static convergence for simply supported plate (b/a=2). 



Next, convergence of the plate natural frequencies is investigated as shown in Fig. 4.4 
for a moderately thick plate (a/h=10). In this study, a full plate model is used so that both 
symmetric and anti symmetric modes are captured. The first five natural frequencies for the 
HOT, which are normalized as indicated in the previous section, are presented for mesh 
sizes ranging from 4x4 to 12x12. These results are compared to the exact values of the 
natural frequencies obtained using FSDT (Reddy, 1984). The 8x8 mesh provides very 
good results and the 12x12 mesh results agree with the exact solutions within 1%. 
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exact (mode 5) 

exact (mode 4) 
exact (mode3) 

exact (mode 2) 
exact (mode 1 ) 


■ 

mode 1 

a 

mode 2 

□ 

mode 3 

B 

mode 4 

§§ 

mode 5 


Mesh size 

Fig. 4.4 Convergence of natural frequencies using HOT for a simply supported plate 

(b/a=l, a/h=10). 


4-JL3 Thick isotropic cantilever plates : The higher order theory is now applied to the 
problem of determining natural frequency parameters for moderately thick and thick 
cantilever plates to determine its accuracy. Dimensions of the cantilever plate are presented 
in Fig 4.5 and the material properties are identical to those used in Section 4.1. The 
frequencies are normalized as follows for comparison with published results. 

.2 


co* = coa' 


(4.1.3) 


where 



(4.1.4) 


First, the higher order theory (HOT) is compared to published results for a moderately 
thick plate (b/h=20). Results using the HOT are presented in Table 4.3 using an 8x8 mesh. 
Results obtained using a Ritz solution with polynomial shape functions and also 
NASTRAN results using a standard triangular element are also presented. Aspect ratios of 
both a/b=l and a/b=3 are shown for the first eight modes of plate vibration. Agreement 



.between the HOT and both the Ritz and NASTRAN results is very good. The largest 
deviation between the HOT and the other approaches occurs for the sixth inode for a/b=l. 
However, the Ritz and NASTRAN results do not agree in this case either so that no 
conclusion can be drawn to the accuracy of any of the methods presented for this mode. 
NASTRAN appears to over predict the frequency for the eighth mode since both the HOT 
and the Ritz approaches agree for this mode. The HOT agrees well with the first seven 
frequencies for the plate with an aspect ratio of a/b=3. This is true even for frequencies 4-6 
which are closely spaced. The last frequency presented (mode eight) deviates slightly from 
the Ritz and NASTRAN solution. Again, the Ritz and NASTRAN results do not agree 
either, so no special significance is placed on this small difference. 

Next, natural frequencies obtained using the HOT (8x8 mesh) for very thick cantilever 
plates are compared to published values for aspect ratios of a/b=l and a/b=2 (Table 4.4). 
These published values are obtained using a three dimensional Ritz approach, again with 
polynomial shape functions, and NASTRAN results using a standard three dimensional 
elements. It must be noted that it can be difficult to choose shape functions for the Ritz 
approach for unusual geometries and the three dimensional NASTRAN results are very 
expensive to obtain due to the large number of degrees of freedom required. The HOT can 
be easily used with complicated geometries with far fewer degrees of freedom since it is a 
two dimensional model. The HOT agrees very well with the Ritz and NASTRAN results 
for both aspect ratios. The largest deviation occurs in the seventh and eighth modes of the 
plate with an aspect ratio of a/b=2. These frequencies are slightly over predicted by the 
HOT. 

The results presented in this section correlate cantilever isotropic plate natural 
frequencies obtained using the HOT and published values. These frequencies agree very 
well in most cases. Small deviations occur in a few cases, but often the published values 
disagree as well. Therefore, confidence is obtained that the HOT correctly predicts 



frequencies for cantilever plates, including very thick plates where shear deformation is 
significant. 



Fig. 4.5 Isotropic cantilever plate. 


Table 4.3 Normalized natural frequency parameters (CD*) for moderately thick isotropic 


cantilever plates (b/h=20) 


a/b 

Mode no. 

HOT 3 

Ritz b 

NASTRAN 0 


1 

3.465 

3.452 

3.475 


2 

8.417 

8.504 

8.547 


3 

21.06 

21.09 

21.42 

1 

4 

26.73 

27.34 

27.27 


5 

30.32 

30.07 

31.24 


6 

43.93 

46.59 

44.42 


7 

52.17 

53.52 

- 


8 

59.59 

59.01 

64.66 


1 

3.233 

3.420 

3.419 


2 

20.95 

21.25 

21.12 


3 

21.33 

22.09 

21.38 

3 

4 

59.99 

60.37 

60.19 


5 

64.11 

65.40 

63.47 


6 

65.87 

67.40 

66.42 


7 

1 18.01 

122.00 

120.40 


(a) 8x8 mesh using current HOT ” ™ 

(b) Ritz solution using polynomial shape functions (Zienkiewicz, 1977) 

(c) NASTRAN using CTRIA2 elements (Ramamurti and Kielb, 1984) 
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Table 4.4 Normalized natural frequency parameters (co*) for thick cantilever plates (b/h=5) 


a/b 

Mode no. 

HOT 3 

Ritz b 

NASTRAN C 


1 

3.3565 

3.3687 

3.3624 


2 

7.4805 

7.3397 

7.3941 


3 

10.985 

10.984 

10.943 

1 

4 

17.890 

17.694 

17.672 


5 

22.644 

23.688 

22.148 


6 

24.47 

24.999 

23.950 


7 

26.185 

26.234 

26.227 


8 

29.733 

29.388 

29.187 


1 

3.384 

3.3397 

3.411 


2 

13.493 

12.459 

13.283 


3 

14.602 

14.390 

14.452 

2 

4 

20.375 

19.596 

20.364 


5 

42.735 

38.771 

41.908 


6 

52.295 

52.282 

52.334 


7 

53.937 

52.469 

53.483 


8 

56.527 

54.814 

54.737 


(a) 8x8 mesh using current higher order theory 

(b) Ritz solution using polynomial shape functions (McGee and Leissa; 1991) 

(c) 14x14x3 mesh using MSC/NASTRAN CHEXA elements (McGee and Leissa; 1991) 
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4.2 Fundamental Correlation of Orthotropic Laminates 

In this section, the higher order theory (HOT) is compared to several existing theories 
to determine the accuracy of displacements, stresses and natural frequencies. Static 
displacements and stresses are correlated first. Next, the accuracy of natural frequencies 
predicted by the HOT are investigated and mode shapes are discussed. 

4.2.1 Static evaluation : The test article for static evaluation is a simply supported plate with 
dimensions shown in Fig. 4.1. However, now a sinusoidal distributed load is applied and 
the plate is constructed of a three ply [079070°] orthotropic laminate with ply thickness t = 
h/3. The aspect ratio of the laminate is a/b=3 and material properties are as follows. 

Ej = 25 GPa, E 2 = 1 GPa, v J2 = 0.25, G 23 = 0.2 GPa, G 13 = 0.5 GPa, G ]2 = 0.5 GPa, 

The static results due to an applied sinusoidal load using the HOT are compared to four 
existing theories for analyzing composite laminates. The first of these theories is an exact 
elasticity solution for cross ply laminates formulated by Pagano (1970) (Exact). The 
second theory is an exact solution to an alternative form of the higher order theory (Reddy, 
1984) (HSDT) . The last two theories are the first order shear deformation theory (FSDT) 
and classical laminate theory (CLT). The transverse displacement predicted by each theory 

is evaluated at the center of the plate and is normalized as follows 

h 3 F 

w* = w(a / 2,b / 2,h / 2)— =f- x 100 

9o a 

The stresses are evaluated and normalized as follows. 

h 2 

Oj * = cjj (a / 2, b / 2, h / 2) — — y 

9o a 


(4.2.1) 



a 2 * = a 2 (a/2,b/2,h/6) T 

qo a 

h 2 

g 6 * = g 6 ( 0,0,h/2)— y 

qo a 


g 4 * = G 4 (a/ 2,0,0) — «■ 

qo a 


h 

g 5 * = g 5 ( 0,b/2,0)— y 

qo a 
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(4.2.2 a-e) 


Results are presented in Figs 4.6 (a-f) for the current approach as well as four existing 
theories. The normalized displacements for all theories agree very well for a thin laminate 
(a/h=100) as shown in Fig 4.6 (a). However, the displacements predicted by the CLT 
depart radically from the true displacement as indicated by the Exact solution for very thick 
laminates (a/h=4). This is due to transverse shear deformation which is ignored in the 
CLT, but is important to consider in thicker laminates. The HOT agrees very closely with 
the HSDT for all thicknesses and differs by a maximum of only 6% from the Exact solution 
for the very thick laminate while the CLT solution under predicts the displacement by 82% 
in this case. The HOT provides good estimates for the inplane stresses, Gj, G 2 , and G6 as 
shown in Figs. 4.6 (b-d). Again, all approaches agree well for the thin laminate (a/h=100) 
although significant differences are observed for the thicker laminates. Surprisingly, the 
well known FSDT noticeably under predicts these stresses for all thicknesses other than 
a/h=100, 44% in the worse case. The HOT agrees very well with the exact solution for the 
alternate form of the higher order theory (HSDT) and is much closer to the Exact solution 
than both the FSDT and the CLT for all inplane stresses. Results for the transverse shear 
stresses are presented in Figs 4.6 (e-f) for the HOT, Exact and FSDT solutions. These 
stresses are assumed to be zero for CLT and are not presented. All solutions presented for 
g 4 in Fig. 4.6 (e) predict similar values at the location where this stress is evaluated. 
However, it must be noted that the FSDT predicts a constant value of this stress through 
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the thickness^and, does, not satisfy the requirement that this stress vanish on the free 

surfaces. In Fig 4.6 (f), the HOT under predicts G 5 compared to the elasticity solution. It 
agrees very well with the HSDT and errors compared to the Exact solution are 
approximately one half of those using the FSDT. The deviation in this case is possibly due 
to the fact that the HOT is a two dimensional theory which is used to represent a three 
dimensional state of stress while the Exact approach is a pure three dimensional approach. 
Over all, the HOT predicts stresses much more accurately than the FSDT or CLT 
approaches compared to the Exact approaches. Very close agreement with HSDT is also 
observed from which it can be concluded that the HOT can be used to accurately predict 
displacements and stresses in more complex laminates where exact solutions are not 
practical. 





4.2.2 Dynamic evaluation : Natural frequencies predicted by the higher order theory (HOT) 
of an orthotropic composite laminate are investigated for a variety of geometries and 
stacking sequences. The test article consists of a four ply orthotropic laminate with 
identical material properties as those used in the previous example. The geometry of this 
laminate is shown in Fig. 4.5 with one end fixed. The stacking sequence of the four ply 
laminate is [±0] s which is defined by a singe variable 0. Each ply is assumed to have 
thickness t=h/4. Normalized frequencies using the HOT are compared to frequencies 
predicted by NASTRAN using two dimensional, shear deformable CQUAD4 elements 
which are based on first order shear deformation theory. Results are also obtained for 
classical laminate theory (CLT) which is obtained by setting the higher order terms in the 
HOT to zero. All of the results in this section, including CLT, HOT and CQUAD4 were 
obtained by the author. 

Results to assess the accuracy of natural frequencies for orthotropic laminates using the 
HOT are presented in Tables 4.5 - 4.8. Frequency parameters for two different aspect 
ratios (a/b=l,2) are presented for several different stacking sequences (0=0°, 15°, 30°, 45°, 
60°) and four different laminate thicknesses ranging from thin to very thick 
(a/h= 100,25, 10,5). Five modes are obtained in each case where the frequencies are 
normalized as follows 


GOj* = coja 2 x 10 


i = 1,2, 3, 4, 5 


(4.2.3) 


where 
D = 


E,lv 


2(l D 12 ^2l) 


(4.2.4) 


Correlation between the HOT and the NASTRAN (CQUAD4) results is very good. In 
general, the HOT predicts slightly higher frequencies compared to CQUAD4. The 
formulation of the CQUAD4 elements requires shear correction factors and reduced 
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integration techniques which are tricks needed to avoid shear locking behavior for ,thin 

elements. These tricks could result in an overly “soft” element for highly orthotropic 
materials, such as used in this example. The mode shapes for 0=0° and a/b=l indicate that 
the first five frequencies correspond to the first bending, first torsion, first camber, second 
camber and second bending modes, respectively. Increasing 0 from 0=0° to 0=30° results 
in a decreased bending stiffness and an increased torsional stiffness. Both the HOT and the 
CQUAD4 results display this trend for the bending and torsional frequencies as shown in 
Tables 4.5 - 4.8. Increasing 0 also results in a more pronounced coupling between the 
bending and torsional behavior of the laminate. This is especially evident in the higher 
modes. Beyond 0=30°, the bending stiffness continues to decrease while the torsional 
stiffness also decreases as indicated by a reduction in the second frequency. Although in 
the static case the maximum torsional stiffness occurs at 0=45°, this is not exactly true in 
the dynamic case since the mode shapes exhibit bending-twist coupling. Although in the 
static case the laminate becomes more stiff in torsion, the influence of the softer bending 
behavior reduces, rather than increases the natural frequencies for these modes. These 
features are even more evident for the case of a/b = 2. Here, the first five frequencies 
correspond to the first bending, first torsion, first camber, second bending and second 
torsional mode shapes, respectively. Again, the distinction between the bending and 
torsional mode shapes becomes blurred for the higher modes as 0 is increased due to the 
coupling of these features. In this case, the frequency corresponding to the second bending 
mode decreases as 0 is increased as expected while the frequency corresponding to the 
second torsional mode increases up to 0=30°, then decreases. According to the classical 
laminate theory (CLT), the normalized natural frequency parameters remain constant for 
increasing laminate thickness since transverse shear deformation is neglected. This leads to 
serious errors for thicker laminates, especially for highly orthotropic materials, as observed 
in Tables 4.7 and 4.8. In the worst case, the CLT over predicts the natural frequency by as 
much as a factor of four (Table 4.8) for a/b=2, a/h=5 and 0=45°. It can therefore be 
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concluded that the CLT is unacceptable tor phsdictlbrf of natural' frequencies in thicker 

orthotropic laminates. However, the HOT agrees well with the NASTRAN results and is 
considered effective for dynamic analysis of composites, even thick and highly orthotropic 


laminates. 
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Table 4.5 Normalized frequency parameters (co*) for cantilever orthotropic plates, a/h=100 


a/b 

0 

Solution 

mode 1 

mode 2 

mode 3 

mode 4 

mode 5 



HOT 

5.580 

6.119 

10.300 

21.567 

34.578 


0 

CQUAD4 

5.544 

5.994 

9.758 

19.799 

33.731 



CLT 

5.590 

6.132 

10.318 

21.605 

34.995 



HOT 

4.725 

7.030 

12.992 

24.528 

30.072 


15 

CQUAD4 

4.660 

6.917 

12.341 

22.388 

28.228 



CLT 

4.734 

7.050 

13.045 

24.631 

30.426 



HOT 

3.537 

8.214 

16.427 

22.825 

28.327 

1 

30 

CQUAD4 

3.500 

8.041 

15.323 

21.544 

26.118 



CLT 

3.546 

8.254 

16.568 

23.043 

28.612 



HOT 

2.281 

8.296 

12.959 

20.964 

27.222 


45 

CQUAD4 

2.261 

8.061 

12.169 

20.030 

25.718 



CLT 

2.289 

8.346 

13.068 

21.120 

27.495 



HOT 

1.404 

6.763 

8.387 

19.655 

23.176 


60 

CQUAD4 

1.400 

6.501 

8.118 

18.642 

22.098 



CLT 

1.407 

6.796 

8.411 

19.767 

23.299 



HOT 

5.579 

7.492 

30.459 

34.662 

37.047 


0 

CQUAD4 

5.543 

7.324 

28.929 

33.786 

35.575 



CLT 

5.589 

7.514 

30.519 

35.079 

37.494 



HOT 

4.537 

10.848 

29.256 

33.452 

41.438 


15 

CQUAD4 

4.489 

10.628 

28.033 

31.873 

39.805 



CLT 

4.549 

10.925 

29.612 

33.658 

41.984 



HOT 

3.201 

13.993 

19.961 

40.237 

46.442 

2 

30 

CQUAD4 

3.148 

13.614 

18.956 

38.408 

44.085 



CLT 

3.221 

14.163 

20.257 

40.746 

47.269 



HOT 

1.982 

11.585 

14.790 

34.888 

43.521 


45 

CQUAD4 

1.910 

11.173 

14.002 

33.485 

40.696 



CLT 

2.003 

11.715 

15.055 

35.450 

44.341 



HOT 

1.115 

7.845 

12.136 

23.140 

36.887 


60 

CQUAD4 

1.260 

7.657 

11.383 

22.322 

34.104 



CLT 

1.123 

7.884 

12.299 

23.326 

37.427 
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Table 4.6 Normalized frequency parameters (co*) for cantilever orthotropic plates, a/h=25 


a/b 

0 

Solution 

mode 1 

mode 2 

mode 3 

mode 4 

mode 5 



HOT 

5.437 

5.937 

10.027 

20.899 

22.926 


0 

CQUAD4 

5.391 

5.805 

9.512 

19.295 

22.849 



CLT 

5.590 

6.132 

10.338 

21.605 

34.995 



HOT 

4.595 

6.755 

12.304 

23.126 

26.003 


15 

CQUAD4 

4.519 

6.630 

11.658 

21.143 

24.188 



CLT 

4.734 

7.050 

13.045 

24.631 

30.426 



HOT 

3.433 

7.751 

14.913 

20.345 

25.236 

1 

30 

CQUAD4 

3.375 

7.551 

13.810 

19.096 

23.363 



CLT 

3.546 

8.254 

16.568 

23.043 

28.612 



HOT 

2.206 

7.773 

11.862 

19.228 

24.352 


45 

CQUAD4 

2.163 

7.500 

11.099 

18.336 

22.981 



CLT 

2.289 

8.346 

13.068 

21.120 

27.495 



HOT 

1.380 

6.411 

8.102 

18.345 

21.689 


60 

CQUAD4 

1.360 

6.115 

7.837 

17.360 

20.635 



CLT 

1.407 

6.796 

8.411 

19.767 

23.299 



HOT 

5.434 

7.170 

21.766 

29.223 

29.884 


0 

CQUAD4 

5.390 

7.010 

21.686 

28.073 

28.790 



CLT 

5.589 

7.514 

30.519 

35.079 

37.494 



HOT 

4.392 

10.051 

25.348 

30.881 

35.741 


15 

CQUAD4 

4.342 

9.844 

24.068 

29.691 

34.151 



CLT 

4.549 

10.925 

29.612 

33.658 

41.984 



HOT 

3.049 

12.629 

17.577 

33.532 

35.432 

2 

30 

CQUAD4 

2.985 

12.272 

16.651 

32.207 

34.027 



CLT 

3.221 

14.163 

20.257 

40.746 

47.269 



HOT 

1.861 

10.681 

13.024 

19.827 

30.453 


45 

CQUAD4 

1.772 

10.314 

12.388 

17.344 

29.082 



CLT 

2.003 

11.715 

15.055 

35.450 

44.341 



HOT 

1.159 

7.544 

10.982 

15.378 

21.530 


60 

CQUAD4 

1.21 1 

7.359 

10.287 

13.337 

20.817 



CLT 

1.123 

7.884 

12.299 

23.326 

37.427 
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Table 4.7 Normalized frequency parameters (co*) for cantilever orthotropic plates, a/h=10 


a/b 

0 

Solution 

mode 1 

mode 2 

mode 3 

mode 4 

mode 5 



HOT 

4.783 

5.148 

8.962 

9.171 

19.567 


0 

CQUAD4 

4.709 

5.010 

8.609 

9.140 

17.326 



CLT 

5.590 

6.132 

10.318 

21.605 

34.995 



HOT 

4.067 

5.769 

10.279 

17.495 

18.764 


15 

CQUAD4 

3.980 

5.650 

9.829 

16.280 

17.577 



CLT 

4.734 

7.050 

13.045 

24.631 

30.426 



HOT 

3.069 

6.437 

11.552 

14.830 

18.316 

1 

30 

CQUAD4 

3.014 

6.284 

10.850 

14.129 

17.283 



CLT 

3.546 

8.254 

16.568 

23.043 

28.612 



HOT 

1.990 

6.445 

9.443 

13.553 

14.939 


45 

CQUAD4 

1.959 

6.243 

8.977 

12.680 

14.418 



CLT 

2.289 

8.346 

13.068 

21.120 

27.495 



HOT 

1.293 

5.506 

7.217 

10.001 

14.791 


60 

CQUAD4 

1.287 

5.304 

7.037 

9.474 

14.266 



CLT 

1.407 

6.796 

8.411 

19.767 

23.299 



HOT 

4.782 

5.962 

8.706 

19.586 

21.075 


0 

CQUAD4 

4.708 

5.849 

8.674 

18.173 

19.667 



CLT 

5.589 

7.514 

30.519 

35.079 

37.494 



HOT 

3.889 

7.958 

16.524 

17.223 

22.878 


15 

CQUAD4 

3.824 

7.836 

16.187 

16.408 

21.939 



CLT 

4.549 

10.925 

29.612 

33.658 

41.984 



HOT 

2.720 

9.644 

13.021 

13.412 

25.144 

2 

30 

CQUAD4 

2.672 

9.480 

12.416 

12.883 

24.347 



CLT 

3.221 

14.163 

20.257 

40.746 

47.269 



HOT 

1.595 

7.933 

8.804 

9.735 

21.952 


45 

CQUAD4 

1.598 

6.938 

8.568 

9.532 

21.056 



CLT 

2.003 

1 1.715 

15.055 

35.450 

44.341 



HOT 

0.976 

6.154 

6.736 

8.540 

17.189 


60 

CQUAD4 

1.156 

5.335 

6.665 

8.312 

16.995 



CLT 

1.123 

7.884 

12.299 

23.326 

37.427 


/ 




Table 4.8 Normalized frequency parameters (o)*) for cantilever orthotropic plates, a/h=5 


a/b 

0 

Solution 

mode 1 

mode 2 

mode 3 

mode 4 

mode 5 



HOT 

3.617 

3.818 

4.585 

7.284 

12.393 


0 

CQUAD4 

3.442 

3.637 

4.570 

7.203 

10.893 



CLT 

5.590 

6.132 

10.318 

21.605 

34.995 



HOT 

3.113 

4.271 

7.863 

10.149 

11.174 


15 

CQUAD4 

3.010 

4.136 

7.722 

10.095 

10.235 



CLT 

4.734 

7.050 

13.045 

24.631 

30.426 



HOT 

2.484 

4.758 

8.205 

10.008 

11.000 

1 

30 

CQUAD4 

2.423 

4.636 

7.850 

9.579 

10.688 



CLT 

3.546 

8.254 

16.568 

23.043 

28.612 



HOT 

1.704 

4.807 

6.776 

6.777 

11.731 


45 

CQUAD4 

1.685 

4.694 

6.340 

6.536 

10.189 



CLT 

2.289 

8.346 

13.068 

21.120 

27.495 



HOT 

1.210 

4.285 

5.001 

5.742 

8.994 


60 

CQUAD4 

1.190 

4.194 

4.737 

5.640 

8.857 



CLT 

1.407 

6.796 

8.411 

19.767 

23.299 



HOT 

3.619 

4.195 

4.353 

12.406 

13.355 


0 

CQUAD4 

3.441 

4.186 

4.337 

10.892 

12.432 



CLT 

5.589 

7.514 

30.519 

35.079 

37.494 



HOT 

3.012 

5.476 

8.262 

11.087 

14.882 


15 

CQUAD4 

2.936 

5.510 

8.204 

10.178 

14.175 



CLT 

4.549 

10.925 

29.612 

33.658 

41.984 



HOT 

2.246 

6.442 

6.706 

8.887 

16.098 

2 

30 

CQUAD4 

2.239 

6.441 

6.526 

8.41 1 

15.451 



CLT 

3.221 

14.163 

20.257 

40.746 

47.269 



HOT 

1.462 

3.966 

6.373 

6.549 

10.976 


45 

CQUAD4 

1.435 

3.469 

6.390 

6.518 

10.668 



CLT 

2.003 

11.715 

15.055 

35.450 

44.341 - 



HOT 

1.057 

3.075 

5.469 

5.748 

8.594 


60 

CQUAD4 

1.097 

2.667 

5.458 

5.911 

8.497 



CLT 

1.123 

7.884 

12.299 

23.326 

37.427 


4.3 Correlation of Piezoelectric Actuation of Beams and Plates 

4.3.1 Correlation with isotropic piezoelectric beam : The higher order theory is applied to 
the problem of determining static displacements of a beam constructed entirely of 
piezoelectric material. The beam is constructed from two plies of PVDF polymeric 
piezoelectric material. Beam dimensions and poling directions of the piezoelectric plies are 
shown in Fig 4.7. Note that the top layer is polarized in the direction of the applied voltage 
and the bottom layer is polarized in the direction opposite the applied voltage. An electric 
field due to the voltage source is applied across the beam such that the top layer expands 
while the bottom layer contracts causing bending moment and a transverse deflection. The 
material properties for this beam-like piezoelectric structure can be assumed to be isotropic 
with values given as follows 

E—2.0 GPa, v—0.30, dsi=23pm/V 

The transverse tip deflection (w) for the piezoelectric cantilever beam is calculated over 
a wide range of applied electric fields which corresponds to a range of applied voltage 0 < 
V < 500 Volts. Results using the HOT are compared with experimental results obtained by 
Lee and Moon (1989) as shown in Fig 4.8. The tip deflections predicted using the HOT 
are in excellent agreement with the experimental results. It must be noted that PVDF has 
the advantage over piezoceramic materials, such as PZT, that significantly higher electric 
fields can be applied without damage to the piezoelectric properties. The results is that 
larger deflections are possible compared to PZT. However, PVDF is much more compliant 
than stiff PZTs. Therefore, the force developed using PVDF materials is significantly 
lower than that possible using piezoceramics. 
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Fig. 4.8 Tip deflection piezoelectric cantilever beam. 


4.3.2 Correlation of piezoelectric actuation of isotropic and orthotropic plates : The higher 
order theory (HOT) is now used to predict displacements and rotations of isotropic and 
orthotropic laminates which are deformed by energizing surface bonded piezoelectric 
actuators. The results using the HOT are compared to experimental results obtained by 
Crawley and Lazarus (1991). The test case consists of a cantilever laminate with 
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dimensions shown in Fig. 4.9. Fifteen pairs of piezoceramic PZT actuators are surface 

bonded to either an aluminum plate or a graphite/epoxy plate of the same dimensions with 
stacking sequences [0°/45%45°] s and [3073070°] s . The thicknesses of the composite 
laminate and the surface bonded piezoelectric actuators are 0.83 mm and 0.25 mm, 
respectively. The bimorph piezoelectric actuator pairs are energized with opposite electric 
fields to produce plate bending. Material properties are presented in Table 4.9 and will be 
discussed in further detail in a moment. 

The HOT finite element model comprises 77 elements and a total of 616 degrees of 
freedom. Overall, good agreement is observed between the HOT and the experimental 
results for both the aluminum and Gr/Ep plates as shown in Figs. 4.10 - 4.12. The 
nondimensional quantities shown include the plate bending (wi), twist (W 2 ) and camber or 
transverse bending (w 3 ) which are defined as follows. 


Wj = M 2 / W 


(4.3.1) 


w 2 =(M 3 -Mj)/W 


(4.3.2) 


w 3 


M 2 - (M 3 + Mj ) / 2 
__ 


(4.3.3) 


where Mi and M 3 are transverse displacements measured at the outer transverse edges, M 2 
is the transverse displacement measured at the center of the plate and W is the width of the 
plate (Fig 4.9). 

As indicated previously, the material properties for the composite laminate and the 
piezoelectric actuators is presented in Table 4.9. The properties used are identical to those 
used by Seeley and Chattopadhyay (1996) indicated as [1] in Table 4.9. It must be noted 
that the composite substructure is an orthotropic material (Gr/Ep T300/976) while the 
piezoelectric material (PZT G-l 195) is considered an isotropic material in the context of 
two dimensional laminate theories. The values used in this work are either taken directly 
from the original research paper by Crawley and Lazarus (1991), or a reasonable estimate 
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Is' 'made' fof properties not listed in the original study. For instance, the piezoelectric 
coupling can be represented by either a single constant (d 3 j), or two experimentally 
determined parameters as described by Crawley and Lazarus (1991). Neither of these 
values used in their study are listed in their paper for the cases presented. Therefore, a 
reasonable guess must be made. However, other authors have taken unusual liberties with 
the values for the material parameters in an effort to obtain experimental verification. The 
values used by several different authors for verification studies are presented in Table 4.9 
when available. Relevant data which is not given in these references is also indicated. As 
an example, both Chandrashekhara and Agarwal (1993) and Ha et al. (1992) use 
significantly different values of the shear moduli (G 23 , G ]3 , and G ]2 ) for their research. 
Only Gi 2 is given in the original work which is different from both of the values used by 
these authors. Detwiler (1995) does not specify the material properties used, but the 
applied electric fields (E 3 ) for the two cases indicated as specified in their work are 
completely different from the original work of Crawley and Lazarus (1991). All of the 
authors indicated in Table 4.9 present results which correlate with the experimental data. 
However, deviations are also apparent. From this discussion, it can be concluded that the 
HOT agrees with the published experimental data as well as can be expected. Deviations 
are the result of imprecise knowledge of the experimental parameters. 



piezoelectric actuators 


W = 
15.2 cm 


L - 29.2 cm 


Fig. 4.9 Cantilever plate with piezoelectric actuators. 
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Fig. 4.10 Nondimensional static mode shapes for aluminum plate. 
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Fig. 4.1 1 Nondimensional static mode shapes for Gr/Ep [0°/45°/-45°] s plate. 
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Fig. 4.12 Nondimensional static mode shapes for Gr/Ep [30°/3070°] s plate. 


Table 4.9 Material properties for 

cantilever plate experimental validation. 


Material 

Authors 



[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

T300/976 Gr/Ep 








Ei (GPa) 


143 

143 

144.23 

150 

* 

150 

E 2 (GPa) 


9.7 

9.7 

9.65 

9 

* 

9 

V12 


0.3 

0.3 

0.3 

0.3 

* 

0.3 

Gi 2 — G] 3 (GPa) 


6.0 

6.0 

4.14 

7.10 

* 

7.10 

G 23 (GPa) 


2.5 

- 

3.45 

2.50 

* 

* 

PZTG-1195 








E (GPa) 


63.0 

63.0 

63.0 

63.0 

* 

63.0 

V 


0.3 

0.3 

0.28 

0.3 

* 

0.3 

d 3 i (pm/V) 


253 

* 

254 

254 

* 

264 

E 3 (Al) (V/mm) 


630 

630 

- 

- 

- 

- 

E 3 ([0°/±45°] s ) (V/mm) 


394 

394 

394 

394 

630 

400 

E 3 ([0730° 2 ] s ) (V/mm) 


472 

472 

- 

- 

755 

480 


[1] Seeley and Chattopadhyay (1996) 

[2] Crawley and Lazarus (1991) 

[3] Chandrashekhara and Agarwal (1993) 

[4] Haetal. (1992) 

[5] Detwiler et al. (1995) 

[6] Koconis et al. (1994) 

(-) nonessential parameter absent 
(*) essentia] parameter absent 



4.3.3 Correlation of composite beam natural frequencies : The current higher order theory 
(HOT) is now applied to the problem of determining the natural frequencies of beam with 
bonded piezoelectric actuators. Two examples are presented to correlate the HOT with both 
experimental data and other analytical techniques. The dimensions of the two beams used 
in this study which have varying boundary conditions are shown in Figs 13 (a and b) along 
with the actuator locations. The material properties for the aluminum beam sub structure 
and the PZT piezoelectric actuators used for both examples are assumed to be isotropic and 
are given as follows 

Al: E-70 GPa v=0.3 p=2700 kg/m 3 PZT: E=63 GPa v=0.3 p = 7600 kg/ m 3 

ds)=254pm/V 

In the first example the first five natural frequencies of a simply supported aluminum 
beam with dimensions shown in Fig 4.13 (a) are determined using the HOT and compared 
to experimental results obtained by Richard and Cudney (1993). A piezoelectric patch of 
dimensions 38x38mm and thickness 0.3 mm is bonded to one surface of the beam at the 
location indicated in Fig 4.13 (a). Although the stiffness of the piezoelectric patch is 
similar to that of aluminum, it is almost three times as dense. Therefore, the natural 
frequencies of the beam with the piezoelectric patch are different from those of a similar 
beam with no piezoelectric patch. Results obtained using the HOT and from the 
experiments are presented in Fig. 4.14 where it is observed that the first five natural 
frequencies of the beam agree fairly well with the experimental results. In general, the 
frequencies are over predicted by the HOT indicating that the structure is slightly softer than 
indicated by the material and geometric properties alone. This is most likely due to 
nonideal boundary conditions which allow small displacements that cause the experimental 
frequencies to be lower than if perfect boundary conditions were possible. The largest 
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deviation occurs for the fifth mode and is less than 10%. Therefore, correlation between 

the HOT and these experimental results is satisfied. 


(a) 


v 96 mm 


/ ■/ 58 mm 


450 mm 


top view 


38 mm 


(b) 


152 mm 


20 mm 


k 


0.3 mmj ' 


side view 


3.17 mm 


1 5.2 mm 


1.52 mm ^ 


T 


A 

T 


[ | aluminum sub structure 

j||| piezoelectric actuator 

Fig. 4.13 Beams with piezoelectric actuators (a) simply supported (b) cantilever. 
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Mode number 


Fig. 4.14 Correlation of HOT with experimentally determined natural frequencies 

The next example is studied to determine how well the HOT predicts the natural 
frequencies of a cantilever beam with a bonded piezoelectric actuator compared to other 
analytical approaches. The test article is again an aluminum beam with a PZT actuator, 
only now cantilever boundary conditions are used. Dimensions of this beam are shown in 
Fig. 4.13 (b) where it is shown that the PZT actuator covers the entire top surface of the 
beam. The first ten natural frequencies are computed using the HOT. These results are 
compared to the same frequencies determined using the first order shear deformation theory 
(FSDT) and the layerwise theory of Reddy (1991) (LT). It must be noted that 
implementation of the layerwise theory requires discretization through the thickness 
direction as well as the inplane directions which improves the accuracy of the theory, but 
significantly increases the computational effort as well. Results obtained using the FSDT, 
HOT and LT approaches are presented in Table 4.10. In general, frequencies predicted 
using the HOT lie between the values predicted by the FSDT and the LT. The FSDT is 
generally slightly stiffer than either the HOT or the LT due to its crude approximation for 
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the transverse shear stress. The LT can be considered the most accurate of the three 

approaches, but is also the most expensive by a significant margin. The frequency for the 
lowest mode is approximately 2% higher than either the FSDT or the LT. For this study, 
both the FSDT and the LT are implemented using 30 lengthwise elements while only 15 
lengthwise elements are used for the HOT. Further mesh refinement for the HOT would 
most likely reduce this effect. Overall, the trends in the frequencies predicted by the HOT 
correlate well with the FSDT and the LT. 


Table 4. 10 Correlation of natural frequencies (Hz) for cantilever beam. 


Mode Number 

FSDT 

HOT 

LT 

1 

537.6 

548.4 

538.4 

2 

3193.0 

3189.1 

3180.7 

3 

7576.6 

7583.7 

7565.6 

4 

8324.0 

8331.3 

8200.9 

5 

14970.1 

15005.4 

14434.2 

6 

22428.1 

21982.5 

21371.3 

7 

22984.4 

22480.6 

22709.0 

8 

31056.9 

30761.5 

28886.5 

9 

37901.0 

36169.6 

33985.3 

19 

39994.7 

38001.4 

37947.8 


4.4 Correlation of Debonded Orthotropic Laminates 

In this section, the higher order theory (HOT) is used to predict displacements and 
natural frequencies of a debonded composite laminate. The geometry of the test structure is 
shown in Fig 4.15. This laminate consists of four plies, each with thickness t=h/4, and the 
stacking sequence is [0°/90°] s . Debonding is introduced at the tip of the laminate between 
the third and fourth plies shown in the side view of Fig. 4.15. The length of the debonding 
(Ld) is described by the nondimensional debonding parameter (3 (P=L d /a). The static 
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displacements are caused by a transversely applied tip load which is uniformly distributed 

across the width of the laminate, but is applied directly to the fourth ply which becomes 
separated from the rest of the laminate due to debonding. The material properties of the 
laminate are as follows 

E i = 60 GPa, E 2 — 25 GPa, V ]2 = 0.25, G 23 = 4.8 GPa, G/j = 12 GPa, Gj 2 = 12 

GPa, p ~ 1500 kg/m 3 . 

The static and dynamic results obtained using the HOT are compared with results 
acquired using NASTRAN. The HOT analysis, which is two dimensional in nature, 
consists of either a 10x10 or 12x12 mesh, depending on the length of the debonding. The 
nondebonded region requires only one element through the thickness while two elements 
through the thickness are used to model the debonded region. Three dimensional CHEXA 
elements in NASTRAN are used to correlated the HOT results. These elements contain 
eight nodes with three degrees of freedom at each node. A 20x20x4 mesh is used in this 
case. This is the smallest practical mesh for the example presented here due to limitations 
on the aspect ratios of the CHEXA elements. It must be noted that many more degrees of 
freedom are required for the NASTRAN results which results in significantly more 
computation time and memory requirements. To obtain correlation for displacements and 
natural frequencies, several different laminate thicknesses are studied ranging from thin to 
very thick (a/h=100, 25, 10) for different debonding lengths ((3=0, 0.10, 0.25, 0.50). 
Both the HOT and NASTRAN results were obtained by the author. 



z 


Top View 


Side view 



[079079070°] 



Fig. 4.15 Debonded cantilever composite laminate with tip load 


4.4.1 Static correlation of debonded composite laminate : The static deflections of the 
debonded composite laminate due to an applied tip load are evaluated first to judge the 
accuracy of the HOT compared to a three dimensional NASTRAN analysis. The transverse 
deflections are normalized as follows 


w^ = w- 


h 3 E : 

q 0 a 4 


(4.4.1) 


In Figs 4.16 (a-d), w* is presented for the case of a moderately thick laminate (a/h=25) 
for varying debonding lengths ((3=0, 0.10, 0.25, 0.50) for both the HOT and NASTRAN 
results. Obviously, the case of [3=0 represents the limiting case for a composite laminate 
with no debonding. The nondimensionalized deflections are presented along the center of 
the laminated from the clamped end to the tip (w*(x,b/2)). Deflections both above and 
below the debonding are shown which occurs as the debonded section pulls away from the 
rest of the laminate. This is because the tip load is applied to debonded ply, which is the 
fourth and outermost ply. It can be seen that the deflections in the regions above and below 
the debonding become largely independent of one another as the debonding length 
increases. Increasing the debonding length also significantly increases the 
nondimensionalized deflections since the structure is weakened by the presence of 
debonding. It must be noted that the deflections of the debonded composite laminate 
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compared to the geometric dimensions of the laminated appear greatly exaggerated in Figs 
4. 16 (a-d) to facilitate determining the accuracy of the HOT. Agreement between the HOT 
and NASTRAN is excellent, even for large debonding lengths ((3=0.25, 0.50) . The HOT 
predicts deflections which are very similar to those predicted by NASTRAN for this 
moderately thick (a/h=25) laminate. 

The HOT is used to predict deflections of debonded composite laminates for several 
different thicknesses ranging from thin to very thick (a/h= 100,25, 10). The tip deflections 
for these three cases are presented in Figs 4.17 (a-c) where the deflections are normalized 
as before. Both the deflection of the debonded ply above the debonding at the tip of the 
laminate and the other region below the debonding are shown for values of (3=0, 0.10, 
0.25, 0.50. Results are presented for these deflections using both the HOT and 
NASTRAN. As before, agreement between the HOT and NASTRAN is excellent. The tip 
deflection of the debonded region increases significantly as the debonding length is 
increased. Since this region is the most compliant region of the structure, is deflects the 
most. The deflection of the region below the debonded region actually decreases as the 
debonding length is increased since it is significantly stiffer than the debonded region and 
the force is applied only to the debonded region. The HOT under predicts the deflection of 
the debonded region by a small amount for large debonding lengths ((3=0.50). This is 
possibly due to the following rational. The stress gradients near the interface of the 
debonded and nondebonded region are large. The fine mesh used for the three dimensional 
NASTRAN analysis, which includes discretization through the thickness direction, more 
accurately captures local stress concentrations in this area. The course mesh used by the 
two dimensional HOT produces a some what stiffer structure which results in slightly 
smaller displacements for this extreme case. The difference between the HOT and the 
NASTRAN deflections is about 6% in the worse case, but is generally much less. The 
HOT also uses significantly fewer degrees of freedom than the NASTRAN analysis. 



Therefore, the HOT is a reliable, but much more efficient, method of studying the effects of 
debonding in composite laminates. 



0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 


x* x* 

(C) (d) 

Fig. 4.16 Nondimensional displacements of debonded composite laminated. 
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4.4.2 Natural frequency determination of debonded composite laminate : The H OX is. now 

used to predict the natural frequencies of the debonded composite laminate being studied 
shown in Fig 4.15. To obtain correlation for the natural frequencies, several different 
laminate thicknesses are again studied ranging from thin to very thick (a/h=100, 25, 10) for 
four different debonding lengths ((3=0, 0.10, 0.25, 0.50). The HOT results are compared 
to NASTRAN results using three dimensional CHEXA elements as before. Again, all 
results were obtained by the author. The natural frequencies for the first ten modes are 
presented in Tables 4. 1 3 -4.1 3. The natural frequencies are normalized as follows. 

t ° i * = “ i i = 1 >2,...,10 (4.4.2) 

Correlation between the HOT and the NASTRAN (CQUAD4) results is very good. In 
general, the frequencies are reduced for all modes as the debonding length increases. The 
only exceptions to this trend occurs for the first mode of each thickness case studied from 
(3=0 to (3=0.10. This is because a 12x12 mesh is used for (3=0 and a 10x10 mesh is used 
for (3=0.10. The courser mesh has a slight effect on the natural frequency determination, 
but is on the order of 1 % and does not have a significant effect on any of the other modes. 
The lowest natural frequency is not significantly effected by the presence of the debonding 
until it becomes very large. This is because the fundamental mode for the nondebonded 
laminate is also the first bending mode which is relatively insensitive to localized changes, 
such as debonding, until it becomes very large. The second mode, which is a twisting 
mode, is also relatively insensitive to the debonding until it becomes very large. It does 
appear to be more sensitive than the first mode to changes in the debonding length, 
however. The third through sixth modes become increasingly more sensitive to the 
debonding length. A significant drop in these frequencies is observed from (3=0.10 to 
(3=0.25 whereas the first two modes remain mostly unchanged until the debonding length 
is increased from (3=0.25 to (3=0.50. This is because these modes, begin to be effected by 
localized changes in the stiffness properties of the structure due to debonding. The higher 
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modes presented, seven through ten, contain many local inflection points and are 
significantly altered by the presence of even relatively small debondings. These mode 
shapes are altered significantly as (3 is increased and incorporate localized deflections due 
the debonding rather than the global deflections contained in the first few modes. Some of 
the frequencies for modes with large debonding are very closely spaced. It must be noted 
that the debonded structure is not symmetric. Therefore, these frequencies do not 
correspond to repeated symmetric modes, but rather distinct mode shapes due the presence 
of debonding. The closeness of their values is mere coincidence. The frequencies 
corresponding to all modes are also effected by the thickness to length ratio (a/h) of the 
laminates. Increasing the thickness of the laminate results in more deformation due to 
transverse shear effects. The effect is lower natural frequencies as observed in Tables 4.1 1 
-4.13. The HOT captures all of the features of the natural frequencies and mode shapes 
found in the more expensive three dimensional NASTRAN analysis. The results of these 
two approaches correlate very well. Therefore, the HOT can be used as a reliable technique 
to predict natural frequencies of debonded composite laminates. 



Table 4.1 1 Natural frequency parameters for debonded composite laminate (a/h-100) 


mode 

number 

analysis 

type 

M 

debonding length 
(3=0.10 (3=0.25 

(3=0.50 

1 

HOT 

0.2427 

0.2465 

0.2449 

0.2068 


NASTRAN 

0.2436 

0.2436 

0.2424 

0.2057 

2 

HOT 

0.4924 

0.4897 

0.4650 

0.3118 


NASTRAN 

0.4839 

0.4771 

0.4580 

0.2995 

3 

HOT 

1.3908 

1.3687 

1.0501 

0.3120 


NASTRAN 

1.3733 

1.3560 

1.0338 

0.3100 

4 

HOT 

1.5379 

1.5363 

1.1460 

0.5009 


NASTRAN 

1.5401 

1.5352 

1.1128 

0.5107 

5 

HOT 

1.9268 

1.8736 

1.1516 

0.5545 


NASTRAN 

1.9084 

1.8665 

1.1171 

0.5208 

6 

HOT 

3.0349 

2.9024 

1.5212 

0.9745 


NASTRAN 

2.9840 

2.8782 

1.4981 

0.9285 

7 

HOT 

3.3196 

3.1790 

1.5965 

1.0725 


NASTRAN 

3.2837 

3.1869 

1.5720 

1.0694 

8 

HOT 

4.2619 

4.2021 

1.7340 

1.2858 


NASTRAN 

4.2817 

4.1994 

1.6242 

1.2905 

9 

HOT 

4.6268 

4.4877 

1.8610 

1.3950 


NASTRAN 

4.6091 

4.4590 

1.8882 

1.4244 

10 

HOT 

4.9103 

4.6435 

2.3755 

1.6188 


NASTRAN 

4.8430 

4.6071 

2.2228 

1.5634 



Table 4.12 Natural frequency parameters for debonded composite laminate (a/h-25) 


mode 

number 

analysis 

type 

11 

o 

debonding length 
P=0.10 p=0.25 

P=0.50 

1 

HOT 

0.2418 

0.2424 

0.2409 

0.2053 


NASTRAN 

0.2433 

0.2433 

0.2421 

0.2047 

2 

HOT 

0.4856 

0.4770 

0.4532 

0.3101 


NASTRAN 

0.4802 

0.4739 

0.4547 

0.2962 

3 

HOT 

1.3693 

1.3410 

1.0450 

0.3102 


NASTRAN 

1.3599 

1.3425 

1.0073 

0.3058 

4 

HOT 

1.5063 

1.4995 

1.1360 

0.4924 


NASTRAN 

1.5281 

1.5232 

1.0851 

0.5065 

5 

HOT 

1.8767 

1.8214 

1.1389 

0.5528 


NASTRAN 

1.8833 

1.8443 

1.0973 

0.5164 

6 

HOT 

2.9297 

2.8016 

1.4892 

0.9712 


NASTRAN 

2.9239 

2.8223 

1.4756 

0.9232 

7 

HOT 

3.1598 

3.1047 

1.5704 

1.0589 


NASTRAN 

3.2216 

3.1226 

1.5501 

1.0581 

8 

HOT 

3.2483 

3.1379 

1.7190 

1.2648 


NASTRAN 

3.4287 

3.4284 

1.5937 

1.2794 

9 

HOT 

4.0596 

3.9949 

1.8071 

1.3733 


NASTRAN 

4.2060 

4.1202 

1.8580 

1.3990 

10 

HOT 

4.3959 

4.2628 

2.3523 

1.6112 


NASTRAN 

4.5061 

4.3500 

2.1889 

1.5545 



Table 4.13 Natural frequency parameters for debonded composite laminate (a/h= 1 0) 


mode 

number 

analysis 

( yp e 

~co 

II 

o 

debonding length 
p=0. 10 P=0.25 

(3=0.50 

1 

HOT 

0.2375 

0.2378 

0.2365 

0.2033 


NASTRAN 

0.2417 

0.2417 

0.2405 

0.2020 

2 

HOT 

0.4578 

0.4516 

0.4304 

0.3054 


NASTRAN 

0.4659 

0.4608 

0.4420 

0.2880 

3 

HOT 

1.2640 

1.2549 

1.0177 

0.3056 


NASTRAN 

1.2992 

1.2828 

0.9399 

0.2956 

4 

HOT 

1.2778 

1.2580 

1.0896 

0.4700 


NASTRAN 

1.3715 

1.3714 

1.0151 

0.4918 

5 

HOT 

1.3658 

1.3544 

1.1042 

0.5448 


NASTRAN 

1.4675 

1.4629 

1.0403 

0.5044 

6 

HOT 

1.6681 

1.6248 

1.2079 

0.9541 


NASTRAN 

1.7779 

1.7471 

1.3694 

0.9055 

7 

HOT 

2.5372 

2.4515 

1.3471 

0.9946 


NASTRAN 

2.6814 

2.5983 

1.4012 

1.0188 

8 

HOT 

2.9460 

2.8346 

1.4861 

1.0483 


NASTRAN 

2.9444 

2.8533 

1.4737 

1.2307 

9 

HOT 

3.2713 

3.1140 

1.6013 

1.2013 


NASTRAN 

3.2741 

3.2734 

1.5107 

1.2985 

10 

HOT 

3.3203 

3.2350 

1.6485 

1.2773 


NASTRAN 

3.6858 

3.6845 

1.7465 

1.3863 
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4.5 Determination of Penalty Factors 

The determination of the numerical values of the penalty factors, p, and p 2 , which are 
used in conjunction with implementation of the continuity conditions, is an important issue 
which needs to be addressed (Section 3.7). The augmented equations of motion (Eqn. 
3.7.11) are found by minimizing the potential and kinetic energies which include the 
penalty terms to satisfy the continuity conditions (Eqns. 3.7.8 and 3.7.9). Minimizing the 
augmented potential and kinetic energies represents minimizing the actual potential and 
kinetic energy functions while satisfying the constraints introduced by the penalty terms. 
To be effective, magnitude of the penalty terms must be comparable to the magnitude of the 
potential and kinetic energies. This is achieved by setting the penalty terms, p, and p 2 ,- to 
the 1 -norms of the stiffness and mass matrices, respectively. In practice, the natural 
frequencies predicted using this penalty approach are relatively insensitive to the actual 
values of p, and p 2 as demonstrated through the example presented next. 

Consider the case for the four ply cantilever composite laminate used in the previous 
section where b/a=l, a/h=25, (3=0.50. Debonding is introduced at the tip of the laminate 
between the third and fourth plies as shown in Fig. 4.15. The 1-norm of the finite element 
stiffness and mass matrices for this case, as predicted by the higher order theory, are on the 
order of 10 8 and 1 0“ 1 , respectively. Selecting p ,= 1 0 8 and p 2 =10' ] yields very good 
correlation with NASTRAN results as indicated in the previous section. The sensitivity of 
the natural frequencies to changes in the values of these penalty factors is determined by 
varying the values of pi and p 2 over a wide range. 

Results for this investigation are presented in Table 4.14. Here, the natural frequency 
parameters, which are normalized according to Eqn. 4.4.2, are presented for the first mode 
of vibration where the values of p, and p 2 are varied over nine orders of magnitude. In 
Table 4.14, values for the first natural frequency for varying p, are presented by column 
while values for varying p 2 are presented by row. It is observed that the normalized 
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frequency parameter for the first mode changes less than 3.5% over five ...orders of 

magnitude for 10 6 < p) < 10 10 . This frequency parameter remains essentially unchanged 
over seven orders of magnitude for 10' 5 < p 2 < 10 2 . Similar trends are observed for 
frequencies corresponding to the higher modes. Numerical ill-conditioning of the system 
matrices is present for values of p, and p 2 larger than those presented in Table 4.14 relative 
to the 1 -norms of the mass and stiffness matrices. Clearly, this study indicates that the 
accuracy of the results is relatively insensitive to a reasonable selection of the penalty terms 
p, and p 2 . 


Table 4.14 First natural frequencies parameters for Gr/Ep [0790*] s beam, (3=0.50. 


Pi 

10 3 

10 4 

10 5 

10 6 

I 0 7 

10 8 

10 9 

10 i0 

10 !! 

P 2 










KT 5 

0.001 

0.040 

0.133 

0.199 

0.205 

0.205 

0.206 

0.204 

0.187 

10- 4 

0.001 

0.040 

0.133 

0.199 

0.205 

0.205 

0.206 

0.204 

0.187 

10” 3 

0.001 

0.040 

0.133 

0.199 

0.205 

0.205 

0.206 

0.204 

0.187 

icr 2 

0.001 

0.040 

0.133 

0.199 

0.205 

0.205 

0.206 

0.204 

0.187 

10 - ] 

0.001 

0.040 

0.132 

0.199 

0.205 

0.205 

0.206 

0.204 

0.187 

10 ° 

0.003 

0.037 

0.122 

0.199 

0.205 

0.205 

0.206 

0.204 

0.187 

I 0 1 

0.001 

0.024 

0.078 

0.195 

0.205 

0.205 

0.206 

0.204 

0.187 

10 2 

0.000 

0.009 

0.000 

0.000 

0.204 

0.205 

0.206 

0.204 

0.187 

10 3 

0.000 

0.003 

0.010 

0.031 

0.000 

0.205 

0.206 

0.204 

0.187 



5. Results for the Higher Order Theory — 

The developed refined higher order theory provides a framework for the analysis of 
composite laminates of arbitrary thickness with surface bonded/embedded actuators and is 
computationally efficient. Numerical results are presented first to investigate strain and 
stress distributions resulting from piezoelectric actuation which are presented for a 
cantilever composite plate with varying thicknesses. Differences between the developed 
higher order theory and the classical laminate theory are discussed. Static results are 
presented for a composite laminate including debonding. The effect of debonding on the 
open loop mode shape of the composite beams used in the experimental investigation are 
examined. Closed loop mode shapes are also presented. 

5.1 Analysis of Stresses and Strains 

In this section, a detailed investigation of the stresses and strains throughout a 
composite laminate is considered. The test article is a cantilever [0707457-45°] s Gr/Ep 
laminate with embedded PZT actuators with material properties listed in Table 5.1. The 
laminate is 25.0 cm in length and 12.5 cm in width (Fig. 5.1). The outermost plies have 
thicknesses 0.050h, the adjacent plies are 0.200h thick, and the remaining plies are 0.125h 
thick where h is the total laminate thickness indicated in Fig. 5.1. The piezoelectric 
actuators are embedded such that they replace the outermost plies in sections of the laminate 
where they are located. Equal but opposite electric fields of magnitude 500 V/mm are 
applied through the thicknesses of each actuator pair to produce bending which results in a 
tip displacement in the positive z direction (Fig. 5.1). All other loads are assumed to be 
zero. The actuators are located centrally on the laminate and the actuator length and width 
are 15 cm and 5.25 cm, respectively. The finite element model comprises 80 equal-sized 
elements with 630 total degrees of freedom. For comparison, the equivalent model using 
the classical laminate theory consists of 450 degrees of freedom. 
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Table 5. 1 . Material properties... 



Gr/Ep 

PZT 

Ei (GPa) 

143 

63.0 

E 2 (GPa) 

9.70 

63.0 

v 12 

0.300 

0.300 

Gi 2 , Gj 3 (GPa) 

6.00 

24.2 

G 23 (GPa) 

2.50 

24.2 

p (xlO 3 Kg/m 3 ) 

1.39 

7.60 

d 3) (xlO- |2 m/V) 

- 

253 


Top View 


Side View 



Fig. 5.1 Cantilever Gr/Ep [070745 °/-45°] s laminate with piezoelectric actuators ; top and 

side views. 

Strains and stresses are calculated throughout the laminate at the elemental gauss points 
used in the stiffness matrix assembly (4x4 integration). In Figs. 5.2 - 5.3, the axial strain 
and stress distributions are presented, through the thickness, close to the center of section 
P-P as indicated in Fig 5.1. This section is chosen to study the effects of induced strain 




actuation where local stresses are expected to be large, without the influence of global 
boundary conditions such as the fixed end. The transverse shear strain distribution (e ] 3 ) is 
presented in Fig. 5.4 at this same location. Results are presented for laminates of length to 
thickness ratio L/h = 10 and 4 which are considered moderately thick and very thick, 
respectively, using both the classical laminate theory (CLT) and the higher order theory 
(HOT). 

Figures 5.2 (a-b) present the strain distribution through the thickness for the moderately 
thick and very thick laminates, respectively. Although the strain distribution is slightly 
nonlinear for the moderately thick laminate (Fig. 5.2 (a)), the refined theory shows that it is 
significantly nonlinear for the very thick laminate (Fig. 5.2 (a)). The CLT does not 
accurately capture this complexity. The strain is larger near the free edge of the laminate 
due to local deformation from the piezoelectric actuation in the outermost plies, but drops 
off more rapidly away from the actuator near the neutral axis as compared to the CLT. 

The stresses due to the piezoelectric actuation are presented in Figs. 5.3 (a-b) for four 
plies (one half of the laminate thickness). The stresses in the other half of the laminate are 
equal, but opposite, due to the symmetry of the laminate and the actuation mode. In these 
figures, the stresses are normalized, ply by ply, to the largest value predicted by CLT. The 
deviation in the stresses predicted by the HOT and the CLT are small in the outermost ply, 
which is the piezoelectric layer. This is because the effect of the additional strain predicted 
by the HOT is offset by the induced strain due to piezoelectric actuation. In the adjoining 
substrate ply, where the induced strain is zero, the stress is significantly larger than that 
predicted by the CLT. However, the stresses in the plies away from the actuator and near 
the midplane are less than those predicted by the CLT due to the localized nature of the 
stresses resulting from piezoelectric actuation. These results indicate that the contribution 
of the higher order terms from the refined displacement field is significant and is important 
in the prediction of the strains and stresses through the thickness of composite laminates, 
particularly near the source of the induced strain actuation. 



The distribution Of" the transverse shear strain ( c 1 3 ) through the thickness of the 
laminate is presented in Fig. 5.4 for both L/h values and shows significant variation 
through the thickness. The maximum value occurs at the midplane of the plate. As 
expected, the strain increases as the thickness of the plate is increased. It must be noted 
that the transverse strains are assumed to be zero in the classical theory. 




(a) (b) 

Fig. 5.2 Axial strain (a) L/h=10; (b) L/h=4. 



ply 4 ply 3 ply 2 ply 1 


1 12 



0.6 0.7 0.8 0.9 1 1.1 1.2 0.6 0.7 0.8 0.9 1 l.l 1.2 

Normalized axial stress Normalized axial stress 


(a) (b) 

Fig. 5.3 Normalized Axial stress (a) L/h=10; (b) L/h=4. 



Transverse Shear Strain (pm/m) 
Fig. 5.4 Transverse Shear Strain (£ 13 ). 


5.2 Static displacements including debonding 

Numerical results are presented for a cantilever Gr/Ep composite plate with embedded 
piezoelectric actuators and a stacking sequence of [0°/9070°]6 and dimensions L = 25cm, W 
= 10cm, and h = 2.286mm. The ply thickness is 0.127mm. Plies numbered 13-16 are 
replaced with a single PZT piezoelectric layer of thickness 0.381mm which represents a 
commercially available piezoelectric material thickness. Material properties are those of 





Gr/Ep and PZT shown in -Table 5.1. A mesh size of 25x10 elements, resulting in 2002 
degrees of freedom, is used for the nondebonded finite element model. The piezoelectric 
layer is energized with a constant voltage of 150 volts. 

Both longitudinal bending and camber deflections result from the piezoelectric actuation 
for the case with no debonding shown in Fig. 5.5. This type of deformation is produced 
because the piezoelectric layer is offset from the neutral axis of the plate which produces 
bending moments about the X and Y axes. Next, a debonding is introduced at the interface 
between the composite substructure and the piezoelectric layer (plies 12 and 13), located at 
the free end of the plate, through the entire width ( 20 cm < x < 25 cm, 0 cm < y < 10 cm). 
The mesh for the composite substructure, located at the midplane of the nondebonded 
portion of the plate, comprises 20x10 elements. Two additional 5x10 meshes are generated 
above and below the debonding at the midplanes of each region. These meshes are offset 
from the midplane of the nondebonded portion of the plate as shown in Fig. 5.6. It is 
important to note that Figs. 5.5 and 5.6 are not to scale so that the different meshes may be 
easily identified. The piezoelectric layer is again energized with a constant voltage of 150 
volts. The static deflection in this case (Fig. 5.6) is significantly different from the first 
case (Fig. 5.5) where actuation produces bending which results in negative transverse 
displacements. The nondebonded portion of the plate bends downward as before, as does 
the region below the debonding. However, the region above the debonding bends in a 
positive transverse direction. This is because the offset of the piezoelectric layer from the 
local midplane in the debonded region is opposite the offset in the nondebonded region. 
This results in a reversal of the bending moment due to piezoelectric actuation and a 
positive transverse displacement. Therefore, the presence of debonding causes the actuator 
to peel away from the substructure. These results indicate that a thorough understanding of 
the effects due to debonding is critical to designing structures with smart composite 


materials. 
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0 

- 0.0001 

- 0.0002 

- 0.0003 

- 0.0004 

- 0.0005 

- 0.0006 

- 0.0007 



Fig. 5.5 Static deflection due to piezoelectric actuation. 


0.001 27 
0.0005 
0 

- 0.0005 
- 0.001 

Fig. 5.6 Static deflection due to piezoelectric actuation including debonding. 

5.3 Open loop frequencies and mode shapes including debonding 

The effect of debonded piezoelectric transducers on the open loop frequencies and 
mode shape of the composite substructure is investigated in this section. The test 
specimens are identical to those utilized in the experimental investigation corresponding to 
beams 1-4 used in the Chapter 6.. These beams have dimensions 30 x 5.3 x 0.194 cm 
with a stacking sequence of [0°/90°]3 S . Piezoelectric transducers with thickness 0.76 mm 
are bonded near the root of the beam for vibration control. Debonding is introduced at the 
interface of piezoelectric transducer #2 (top surface) and the Gr/Ep substructure at the edge 
of the transducers nearest the root and extends through the width of the beam. As before, 
the parameter (3 represents the length of the debonding. Values of (3 = 0.0, 0.06, 0.12, 
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.0. 18 are investigated for: several different modes. The results obtained using the higher 
order theory (HOT) for the nondebonded cases comprise a 17x3 finite element mesh and 
additional elements are added as needed to account for varying lengths of debonding. 

The open loop mode shapes with and without debonding of piezoelectric transducer #2, 
obtained using the present approach, are presented in Figs. 5.7 and 5.8 and the frequencies 
are shown in Table 5.2. The first five mode shapes in the absence of debonding are shown 
in Fig. 5.7. These mode shapes are typical of the structure being analyzed. The first two 
modes represent bending modes, the third and the fourth modes are twist modes and the 
fifth mode represents the third bending mode as indicated in the first column. The 
frequencies corresponding to these nondebonded ((3=0) modes are presented in Table 2. 
The sequential mode number of each of the frequencies is shown in parenthesis. As the 
debonding length is increased to (3 = 0.06, the magnitudes of the frequencies decrease 
slightly. This is due to a small reduction in the structural stiffness due to debonding. 
However, no significant changes are observed in the mode shapes. 

Further increase in the debonding length produces significant changes in the mode 
shapes. The first several open loop mode shapes for |3 = 0.18 are shown in Fig. 5.8. It 
must be noted that in Fig. 5.8, the wire frame elements used to represent the mode shapes 
are plotted at the midplane of each element. The midplane corresponding to the debonded 
portions of the structure is different from the midplane corresponding to the nondebonded 
structure. Therefore, the elements in the debonded regions appear to be disconnected from 
the rest of the structure. This presentation of the mode shapes is chosen for clarity. Some 
of these modes, such as the first bending mode (Fig. 5.8 (a)) exhibit global deformation of 
the substructure. Other modes, such as the second bending mode (Fig. 5.8 (c), clearly 
indicates local deformation caused by the debonding of the actuators. Therefore, it is 
necessary to study the modes with debonding in greater detail. In Table 2, frequencies 
corresponding to modes which exhibit bending or twisting are indicated using B or T, 
respectively. The letters G and L are used to indicate global or local behaviors, respectively 



and S and A refer to symmetric and anti sy mmetric modes respectively.. JFor. instan;ce,.Figs. 
5.8 (a) and (b) show mode shapes which have characteristics of the first bending mode 
(Bl). In Fig. 5.8 (a), the deformation is primarily global in nature (G). However, Fig 5.8 
(b) indicates that the deformation is predominately local (L) due to the debonding of the 
actuator. The mode shape presented in Fig. 5.8 (a) shows that both the deformation of the 
debonded region and the rest of the substructure are in phase indicating a symmetric (S) 
mode. This is in contrast to the mode shape presented in Fig. 5.8 (b) where the local and 
global deformation are out of phase indicating an antisymmetric (A) mode. Therefore, the 
mode shapes presented in Figs. 5.8 (a) and (b) are denoted GB1S and LB1 A, respectively. 
The frequencies corresponding to these modes are sequential as indicated in Table 2, but 
their values differ significantly (24.208 Hz for LB 2 A and 55.911 Hz for GB2S for 
(3=0.18). 

The presence of debonding can have a significant and counter intuitive influence on the 
values of the frequencies corresponding to the global and local counterparts of modes 
displaying similar physical characteristics. It is observed that increase in the debonding 
length introduces local modes. For example, LB IS ((3=0.12) and LB1A ((3=0. 18) 
correspond to local bending modes which are absent for smaller values of (3 (Table 2). 
Similar observations are made for the higher modes as well. It is also interesting to note 
that although there is a general decrease in the values of the natural frequencies as the 
debonding length is increased, the value of the frequency corresponding to the second 
global bending mode increases (from 1 14.85 Hz to 1 18.12 Hz) as (3 is increased. Another 
point is that for (3=0.12, although the local twist modes are present (LT1A, LT2S), the 
second local bending mode is absent. When the second local bending mode (LB2S) does 
appear at (3=0.18, it is at a significantly higher frequency that either of the twist modes as 
indicated by the magnitude of the frequency and the sequential mode number of the 
frequency in parenthesis The global and local counterparts of the first twist modes for 
(3=0.12 and (3=0.18 display symmetric/antisymmetric behavior while both the global and 



local counterparts of the- second twist modes are symmetric. These observations indicate 
that the dynamic characteristics of the composite laminate can be greatly altered by the 
presence of debonding of the piezoelectric transducers. Therefore, careful attention must 
be paid to the existence of such imperfections in predicting the dynamic response of such 


smart structures. 
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(e) GB3 

Fig. 5.7 (a-e) Mode shapes with no debonding obtained using HOT. 
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(a) GB 1 




(e) GB3 

Fig. 5.8 (a-e) Open loop mode shapes with including debonding obtained using HOT. 



Table 5.2 Change in open loop frequencies due to debonding (Hz)*. 


■$939 1 


(3 = 0.0 

& 

= 0.06 


(3 = 0.12 

h 

= 0.18 

1st 

25.749 

0) 

GB1 

24.958 

(1) 

GB1 

24.555 

(1) 

GB IS 

24.208 

0) 

GB IS 

bending 







128.70 

(3) 

LB IS 

55.91 1 

(2) 

LB1 A 

2nd 

118.07 

(2) 

GB2 

117.01 

(2) 

GB2 

114.85 

(2) 

GB2 

118.12 

(3) 

GB2S 

bending 










382.29 

(9) 

LB 2 A 

1st 

153.28 

(3) 

GT1 

146.26 

(3) 

GT1 

139.42 

(4) 

GT1S 

122.55 

(4) 

LT1S 

twist 







222.56 

(5) 

LT1 A 

146.66 

(5) 

GT1A 

2nd 

311.47 

(4) 

GT2 

299.02 

(4) 

GT2 

266.19 

(6) 

LT2S 

275.84 

(6) 

LT2S 

twist 







302.81 

(7) 

GT2S 

297.55 

(7) 

GT2S 

3rd 

345.54 

(5) 

GB3 

341.98 

(5) 

GB3 

347.48 

(8) 

GB3 

327.82 

(8) 

GB3 


bending 


* (G) global, (L) local, (B) bending, (T) twist, (S) symmetric, (A) antisymmetric 


5.4 Closed loop mode shapes. 

In this section, the effect of the closed loop control due to piezoelectric actuation is 
investigated. The test articles are identical to the nondebonded [0790°]3 S and [457-45°]3 S 
composite beams with bonded piezoelectric transducers (|3=0). The closed loop control 
consists of an accelerometer located at the tip of the cantilever beams and piezoelectric 
transducer #2 (top) is again used as an actuator. The first mode shape for the open and 
closed loop cases of each beam are presented in Figs 5.9 (a-b). In each case, the modes are 
normalized to unity for comparison. The gain for the closed loop case is 8700 which was 
the largest gain considered in the experimental investigation. It is observed from Figs 5.9 
(a) and 5.9 (b) that the closed loop mode shapes are nearly the same as the open loop mode 
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shapes. The control forces act to stiffen the beams near the piezoelectric transducers. 
However, the transducers already add stiffness to the beams in the same locations from 
their material properties and the additional thickness which they provide. Therefore, the 
additional stiffness provided by the closed loop control does not significantly change the 
mode shapes. The higher mode shapes exhibit more localized deformation, such as that 
caused by the piezoelectric actuation. However, the feedback control is most responsive to 
the first mode due to the presence of the low pass filter and less responsive to the higher 
modes. Therefore, the higher closed loop mode shapes do not differ by a significant 
amount either from their open loop counterparts. 


open loop 


open ioop 

• closed loop 


• closed loop 



x x 

(a) (b) 

Fig. 5.9 First open and closed loop mode shapes for (a) [0°/90°]3 S and (b) [45% 

45°]3 s beams. 





Experimental Investigation _ 

The principal objective of the current research is to develop practical analysis tools for 
studying composite structures incorporating piezoelectric sensors and actuators. However, 
it is necessary to validate these tools through experimental correlations. Furthermore, no 
published experimental results have been presented to date which address the issue of 
partially debonded piezoelectric actuators. Many practical issues, which are not obvious in 
the development of the mathematical model, must also be addressed. Therefore, an 
experimental investigation was performed as part of the research for this dissertation. The 
goals of this investigation are twofold: 

1. Assess the accuracy of the higher order theory to model composite structures 
incorporating piezoelectric actuators including the effects of partially debonded 
actuators 

2. Address practical challenges related to the implementation of smart structures. 

In the following sections, details are presented to describe the construction and testing of 
composite structures incorporating piezoelectric sensors and partially debonded actuators to 
reduce vibration. 

6.1 Test Specimen Construction 

The test specimens were constructed to represent variations in: (a) composite ply 
stacking sequence and (b) amount of piezoelectric actuator debonding. Construction 
occurred in two phases. The first phase consisted of the layup and curing of the composite 
specimens. In the second phase, the piezoelectric transducers and accelerometers were 
bonded to the composite substructures. These two phases of construction are described 
next. 

6.1.1 Composite substructure construction : The material selected for the composite 
substructure was HYE-3574 OH Graphite/Epoxy fabric. Two different ply stacking 
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sequences, were studied, fQ790°]3 S and [457-45°]3 S . These ply layups represent commonly 

used laminates in the industry. Two large flat plates (122 x 30 cm), one for each stacking 
sequence, were constructed from the prepreg fabric. The three stage curing cycle was 
standard for this type and geometry of specimen and was specified as follows. 

Stage 1: Heat to 126° C, pressurize to 690 kPa 
Stage 2: Maintain 126° C, 690 kPa for 120 min. 

Stage 3: Cool to 60° C, 690 kPa for 60 min. 

Once the plates were removed from the autoclave, they were cut into strips approximately 
56.26 x 5.34 cm. These strips were of the exact width as the final specimens, but were 
longer than the test articles to allow for clamping at one end. The average thickness of the 
specimens was 1.94 mm. Therefore, the ply thickness was approximately 0.161 mm. The 
average density was calculated by dividing the average mass of each test specimen by its 
volume and was found to be 1507 kg/m 3 . 

6.1.2 Incorporation of piezoelectric transducers .: The second phase of construction 
consisted of attaching the piezoelectric transducers to the composite samples. In this phase, 
it was first necessary to locate practical piezoelectric transducers to be used as sensors and 
actuators. Lead zirconate titanate (PZT) piezoceramic transducers were found to be the 
most useful for the current application of vibration control since they exhibit a relatively 
high elastic modulus and piezoelectric coupling coefficients. Although wafers of raw PZT 
can be purchased, it was impractical to use these in the current investigation for two 
reasons. First, electrodes must be soldered to the material without destroying the 
piezoelectric properties. Second, the bare electrodes would be exposed to the conductive 
Graphite/Epoxy substructure necessitating an insulating layer and thereby complicating 
fabrication. However, PZT transducers were available in a prepackaged format with pre- 
attached leads and electrodes encased in an insulating Kapton package. ACX QP40N 
piezoceramic transducers were selected with dimensions 10.16 x 2.54 x 0.0762 cm. These 



transducers contained two stacks of two PZT wafers shown in Fig. 6.1 and were used as 
both sensors and actuators. 


4 
3 

I 
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Fig. 6.1 ACX QP40N piezoceramic transducer. 

Once practical transducers were obtained, it was then necessary to bond them to the 
composite specimens. Reliable bonding was essential to ensure the most efficient transfer 
of strain from the actuator to the composite substructure and vice versa. Creating a 
bonding layer which was stiff, thin, and free of air pockets was critical to maximize the 
effectiveness of the actuators. The adhesive material was chosen to be Ecobond 45 clear 
epoxy adhesive. It was a stiff, room temperature cure epoxy which had a working life of at 
least an hour and was recommended by the manufacturer of the piezoelectric transducers 
(ACX). 

Several prototypes were constructed to debug the process of bonding the actuators to 
the composite substructure. The composite specimens and piezoelectric transducers were 
quite expensive. Therefore, the prototypes were constructed from plexiglass, which was 
used to represent the composite specimens, and aluminum strips, which represented the 
piezoelectric transducers. Since the plexiglass was clear, the additional benefit of visual 
inspection of the bonding layer from the bottom was possible. The manufacturer of the 
piezoelectric actuators recommended vacuum bagging their products to the host structure 
while the adhesive layer cures. After several attempts, this procedure was determined to be 
impractical since the epoxy was consistently drawn out from the interface of the transducer 
and substructure. A more practical procedure was found to be simply placing weights on 
top of the transducers while the epoxy cured. However, this approach lead to the problem 
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of minuscule air pockets contaminating the bonding layer. Placing the epoxy in a vacuum 

chamber for several minutes before applying it to the bonding surfaces, which allowed the 
air pockets to expand and burst, produced favorable results. Teflon tape was placed on the 
composite substructure around the area to be bonded to prevent excessive flow of epoxy. 
Release paper was placed in between the test specimen and the weights to prevent bonding 
of the weights to the test specimens. Teflon tape, with thickness 5.08 pm, was placed 
underneath portions of the piezoelectric actuators to create pre-existing debonding of the 
actuators. The tape was coated with a thin layer of petroleum jelly to ensure that no 
adhesion between the actuator and composite substructure occurred. After allowing 24 
hours for the epoxy to cure, the specimens were released and the excess epoxy trimmed. 
The entire configuration is presented in Fig. 6.2. Additionally, a single Endevco 2250A 
Micro miniature accelerometer was bonded to the tip of each of the test specimens using a 
cyanoacrylate adhesive. 


x ^ ^ ^ $>; p: ^ ^ ; 
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Fig. 6.2 Bonding configuration for the test specimens. 


6.1.3 Test specimen configuration : Eight composite test specimens were constructed with 
varying stacking sequences and debonding lengths. Piezoelectric transducers were bonded 
to the top and bottom surfaces of the composite substructure as indicated in Fig. 6.3 along 
with a single accelerometer at the tip. Wooden blocks were clamped at one end to provide 





fixed boundary conditions. The test specimens were clamped sideways in a vise during 
testing to further ensure satisfaction of the fixed boundary conditions and minimize 
extraneous vibration from corrupting the data. Debonding was introduced between the top 
piezoelectric transducer and the composite substructure through the width of the transducer 
near the clamped end as shown in Fig. 6.3 (a) The nondimensional debonding length, p, 
is defined as follows. 

P = ^ (6.1.1) 

where L D is the actual debonding length and L is the length of the test specimen. The 
physical dimensions of the test specimens are presented in Table 6.1 and the exact 
configuration of each beam is presented in Table 6.2. Material properties are presented in 
Table 6.3. It must be noted that the d 3 i value listed for the PZT transducers is that of the 
piezoelectric material only. The transducers are really composite structures containing 
several different materials, such as the electrodes and the Kapton casing. Therefore, it is 
more useful to express the piezoelectric coupling in terms of micro-strain per volt for the 
entire composite transducer. This quantity also has a small but noticeable dependency on 
the applied electric field. The value for the transducers used for beams 1-4 ([0790°] 3s ) was 
0.81 m while the value for beams 5-8 ([457-45°] 3s )was 0.73 / m due to the 


lower fields used during testing. 
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Top View (b) 

Fig. 6.3 Test specimen configuration. 


Table 6.1 Test specimen dimensions (cm.). 


L 

Li 

l 2 

W 

t 

h 

30 cm 

3.4 cm 

1 3.7 cm 

5.3 cm 

0.76 mm 

1 .94 mm 


Table 6.2 Test specimen configuration. 


beam 

stacking sequence 

p 

1 

[0790° ] 3s 

0 

2 

[0790°] 3s 

0.06 

3 

[0790°] 3s 

0.12 

4 

[0790°] 3s 

0.18 

5 

[457-45°] 3s 

0 

6 

[457-45°] 3s 

0.06 

7 

[457-45°] 3s 

0.12 

8 

[457-45°] 3s 

0.18 



Table 6.3 Test specimen material properties. 


HYE-3574 OH Gr/Ep 

PZT 

Ei (GPa) 

114 

6.9 

E 2 (GPa) 

9.5 

6.9 

G] 2 =G 13 (GPa) 

4.7 

2.6 

G 23 (GPa) 

2.1 

2.6 

U12 

0.30 

0.31 

p (kg/m 3 ) 

1507 

5000 

d 3 i (pm/V) 

- 

-179 

k 3 

- 

1800 


Wires had to be attached to the piezoelectric transducers and the accelerometer on each 
specimen. This was accomplished via two channels in the clamped wooden blocks (Fig. 
6.4). Digi-Key flex cables plug into the piezoelectric transducers and run through the 
channels in the wooden blocks. The same was true for the wires which were connected to 
the accelerometers. The flex cable was then attached to AWG 22 stranded core wire which 
lead to a terminal bus with screw connectors. This intricate wiring scheme was necessary 
for two reasons. First, it allowed for the least possible interference with the clamped 
boundary conditions. Second, the four pin connectors mounted on the piezoelectric 
transducers were extremely fragile. The extension to the terminal bus provided a rugged 
and convenient connection for practical test conditions. 



Fig. 6.4 Test specimen wire configuration. 
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6.2 System Identification 

The main goal of the experimental investigation is to correlate the developed 
mathematical model with the composite beams which were constructed. This requires 
experimental determination of the frequencies and damping ratios for each test specimen. 
These quantities are extracted from the raw experimental data using a stochastic method 
known as the Auto Regressive Moving Average with eXogenious input model (ARMAX) 
(Lee and Fassois; 1990: Mignolet et al.; 1993: Mignolet and Red-Horse; 1994). Using this 
method, the system is assumed to be linear and represented in the continuous domain by 
the following standard set of second order ordinary differential equations. 

Mx + Cx + Kx = f(t) (6 2 1) 

where M, C and K are the standard mass, damping and stiffness matrices, respectively, x 
is the displacement vector and f(t) is the time varying force vector. This system is 
parameterized in the discrete domain as follows 

A(z)x(z) = B(z)f(z) + C(z)e(z) (6.2.2) 


where the term on the left hand side represents the output, the first term on the right hand 
side represents the input while the second term represents the noise. The parameters 
contained in A, B and C are determined using the ARMAX system identification 
technique and z is the unit lag operator. The system poles, Zj, are determined from the 
solution of the following characteristic equation. 

det(A(zi)) = 0 (6.2.3) 


The natural frequencies and damping ratios of the structure are obtained by mapping the 


system poles from the discrete to the continuous domain as follows. 

i 


CO; 


i — ( ln ( r i) +0 > 2 


X 


(6.2.4) 



£i =-ln(ri)( | n(r i ) 2 + ef) 2 
where 

rj = |zj |, 9j = tan 


\ Re ( z i )J 
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(6.2.5) 

( 6 . 2 . 6 ) 


6.3 Open Loop Structural Response 

63.1 Open loop transient response : The open loop natural frequencies and damping ratios 
of the test specimens are experimentally determined in this section The ARMAX system 
identification technique was used on the transient response records of each of the eight test 
specimens. The natural frequencies and damping ratios corresponding to first two modes 
of vibration, which were the first two bending modes, were then compared with those 
predicted by the developed finite element model implementation of the higher order theory. 
The next two modes, which were twisting modes, were neither observable or controllable 
using the current sensor/actuator configuration. Therefore, they were not considered. The 
third bending mode occurred at a frequency nearly three times higher than the second 
bending mode and was not considered either. 

The test setup consisted of clamping each specimen in a vise which provided fixed 
boundary conditions. The tip of each beam was lightly tapped and the charge output of 
piezoelectric transducer #1 (bottom) was converted to a voltage and recorded using the data 
acquisition system. The entire configuration is shown in Fig. 6.5 



Fig. 6.5 Open loop transient response configuration: 
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The response was measured in terms of a voltage which is directly proportional to the 

charge generated by piezoelectric transducer. This charge is proportional to the strain in the 
piezoelectric transducer and therefore represents the displacement of the beam at a given 
time. The voltage history for beams 1 and 5 ([0°/90°]3 S and [45%45°] 3s , (3=0) is presented 
in Figs 6.6 and 6.7. 



0 0.5 1 1.5 2 2.5 3 3.5 4 

Time (seconds) 

Fig. 6.6 Open loop transient response beam 1 ([0°/90°]3 S , (3=0). 
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Fig. 6.7 Open loop transient response beam 5 ([-45°/45°]3 S , [3=0). 


A fast Fourier transform on the voltage response records of beams 1 and 5 reveals the 
frequency content. These transforms are presented in Figs. 6.8 (a-b) 



(a) (b) 

Fig. 6.8 FFT of open loop transient response (a) beam 1 (b) beam 5 


The first peak in Figs 6.8 (a-b) represents the frequency of the first mode while the second 
peak represents the second mode. From Fig. 6.8 (a), the frequencies of the first two 
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modes for beam 1 were approximately 25 Hz and 120 Hz while the first tw.aJxequencies of 

beam 5 were approximately 18 and 65 (Fig. 6.8 (b)). The open loop natural frequencies of 
beams 1-4 were expected to be similar since all of these beams have the same stacking 
sequence. The same was true for beams 5-8. Some small variation was expected due to 
the debonding, though. Once the frequencies had been approximately determined, 
damping in each of the beams was estimated using the logarithmic decrement approach as 
follows. 



(6.3.1) 


where 6 is the logarithmic decrement and x n is the amplitude after n cycles. The damping 
ratio, is approximated as follows 



Using this approach, the damping ratio for the first mode in each beam, which is the 
dominate mode in the transient response, was found to be quite small (^=0.0025). In 
general, the decay envelope of the structure is expressed as follows 

x(t) = e"“ & (6.3.3) 


Starting from the initial time to, the fraction of amplitude A which remains after some time t 


is as follows. 


Ae -<(u+0 

Ae-^ l ° 


= e - 0 * 1 


(6.3.4) 


In order to accurately measure the damping ratios of the test specimens, a significant decay 
in the vibrational amplitude had to be obtained. The amount of time required to allow for 
the structure to decay to a certain level was found by rearranging the above equation as 


follows. 



v..,. .... .... ... . ■ . . 135 

Engineering judgment suggested that the vibrational ampl i tudc shou Id be reduced 'to less ' 

than 10% of its original value to allow for accurate determination of the damping ratio. 
Using the above equation, along the lowest expect natural frequency (18Hz) and damping 
ratio (0.0025), the time required was approximately eight seconds. Therefore, eight 
seconds of data were recorded to determine the open loop transient response of each beam. 

Once the transient voltage histories of each beam were obtained, the ARMAX technique 
was applied to the voltage records to more accurately determine the open loop natural 
frequencies and damping ratios for the first two modes of each test specimen. Model 
orders ranged from five to twenty five. Lower model orders were consistently inaccurate 
while higher model orders resulted in numerical problems due to over determination. 
Frequency estimates were consistent over the range of model orders considered. However, 
the damping ratio estimates varied somewhat due to the presence of nonlinear damping. 
Therefore, several different models, using different model orders, of each voltage history 
were generated. Then, the average damping ratios were calculated. The sampling rate used 
was 250 Hz. This is more that two times the highest expected frequency to be determined 
for accuracy, but is not so high as to degrade the signal to noise ratio of the data. 

The experimentally determined frequencies (EXP) for the first two modes are presented 
in Table 6.4 and are compared to the same frequencies predicted using the higher order 
theory (HOT). In general, the agreement is very good. Although the mass of the 
accelerometers is very small (0.4 g), it is included in the mathematical model. The 
frequencies were expected to decrease slightly with increasing debonding length which was 
observed as a general trend. Since the change is relatively small, it is not clearly observed 
in the experimental results of the first mode of the [0790°] 3s beams due to experimental 
uncertainty, but it is observed in the [457-45°] 3s beams. However, this trend is always 
observed in the HOT results. An unexpected change occurred in the frequency of the 
second mode as the debonding length increases from (3=0.12 to (3=0.18. Although the 
structure became less stiff due to the increased debonding length, which would indicate a 
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reduction in the frequencies, the experimental value of the frequency of the second mode 
actually increased. This increase is also observed in the HOT results for both beams and in 
the experimental results for the [0790°]3 S beams as well. The reason for this anomaly is 
that a new localized mode was introduced between the first and second bending modes due 
to the flapping motion of the debonded actuator. This new mode alters the dynamic 
response of the beam which results in a slightly higher frequency for the second bending 
mode. 


Table 6.4 Open loop natural frequency estimates (Hz) 


beam 

EXP 

mode 1 
HOT 

% error 

EXP 

mode 2 
HOT 

% error 

1 

25.1 

25.4 

1.0 

120.6 

116.3 

3.6 

2 

24.5 

24.6 

0.3 

118.7 

115.3 

2.9 

3 

24.6 

24.2 

1.7 

119.4 

113.1 

5.3 

4 

24.5 

23.8 

2.7 

120.3 

116.3 

3.3 

5 

17.3 

16.4 

5.3 

65.5 

68.0 

3.9 

6 

16.3 

15.6 

4.1 

66.3 

67.3 

1.5 

7 

15.5 

15.2 

2.1 

66.6 

67.2 

0.9 

8 

15.4 

14.8 

4.1 

65.9 

68.4 

3.8 


The experimentally determined open loop damping ratios are presented in Table 6.5 
The damping in all cases was found to be less than 1% which was quite small. The 
damping increases from the first to the second mode for beams 1-4 ([0790°] 3s ) while it 
decreases for beams 5-8 ([45745°] 3s ). The damping was generally higher for the 
[45°/45°]3 s beams compared to the [0790°] 3s beams. This was because the epoxy matrix, 
which has inherently higher damping characteristics than the graphite fibers, had a more 
significant effect on the structural response of the [45745°]3 S beams. 
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Table 6.5 Open loop damping ratio estimates 


beam 

EXP- mode 1 

EXP - mode 2 

1 

0.0030 

0.0039 

2 

0.0023 

0.0060 

3 

0.0035 

0.0038 

4 

0.0036 

0.0040 

5 

0.0054 

0.0049 

6 

0.0053 

0.0046 

7 

0.0068 

0.0046 

8 

0.0085 

0.0052 


6.3.2 Open loop frequency response : The open loop frequency response of the test 
specimens is determined next. A sinusoidal voltage with varying frequency is applied to 
piezoelectric transducer #1 and the voltage response of piezoelectric transducer #2 was 
recorded. The experimental setup is presented in Fig. 6.9. A voltage controlled oscillator 
circuit, which was designed, tested, and built, was coupled with the data acquisition 
system and provided the frequency sweep. Acquiring the frequency response data was 
very difficult due to the range of amplitudes at different frequencies. The response of the 
specimens near their respective natural frequencies was very large since the structural 
damping was very small. Therefore, the input voltage, V, was required to be very small to 
keep the structural response within reasonable limits at the resonant frequencies. However, 
the response of the structure at other frequencies was extremely small at the same voltage 
and often difficult to distinguish from ambient noise in the structure and the charge 
amplifier. 
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Fig. 6.9 Open loop frequency response experimental setup. 

The frequency response of test specimens 1 ([0°/90 o ] 3s ) and beam 5 ([45%45 0 ] 3s ), 
which have different stacking sequences but no debonding, are presented in Figs 6.10 (a- 
b). The magnitudes in these plots are relative and are only relevant for comparing the 
higher order theory (HOT) to the experimental data (EXP). The frequency response of test 
specimens 2-4 and 6-8 were similar to test specimens 1 and 5, respectively and were 
largely unaffected by the presence of debonded transducers, except for a small shift in the 
resonant frequencies. Therefore, only the response of test specimens 1 and 5 are 
presented. 

The location of the peaks in the structural response observed by piezoelectric transducer 
#2 correlated very well with the peaks which were observed experimentally for all of the 
test specimens. This was expected since the natural frequencies predicted by the higher 
order theory correlated well with the experimental results as indicated in the previous 
section. The first two modes were bending modes and were observable by the sensor. 
The next two modes were twist modes and were not observable or controllable using the 
current sensor and actuator configuration. The next observable mode, which was the third 
bending mode, has a frequency much higher than the first two modes and was not present 
in the frequency range of interest. Due to the close proximity of transducers #1 (actuator) 
and #2 (sensor), the deformation seen by transducer #2 included both the global response 
of the structure as well as the local deformation produced by transducer #1. The phase of 
each these components below resonance was 0° and they became additive. However, each 
these components above resonance was 180° out of phase and they had a canceling effect 
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on each other. This explains the nonsymmetric shape of the peaks at each of the resonant 
frequencies which was predicted by the higher order theory and verified experimentally. 
(Figs. 6. 10 (a-b)). Since a significant portion of the deformation observed by piezoelectric 
transducer #2 was local deformation caused by the actuation of transducer #1, stability 
problems were encountered for the closed loop control which required information about 
the global response of the structure. Therefore, the accelerometer at the tip was used as a 
sensor instead for the closed loop feedback control presented next. 
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6.4 Control Design 

In this section, topics relating to design and implementation of the feedback control 
system are discussed. A simple analysis based on classical control theory is developed to 
understand basic concepts relating to the control system. The full state space equations are 
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then introduced and the corresponding circuit design for implementation of the controls is 
presented. 

6.4,1 Classical approach to control design : Consider a cantilever beam with a piezoelectric 
transducer used as an actuator and an accelerometer used as a sensor. The equation of 
motion for the first bending mode of the beam, which is assumed to be the most important 
mode, is given as follows. 

q + 2^ b co b q + ojgq = 4 > t (f + F p ) (6.4.1) 

where q, <j), and cob are the modal participation factor, mode shape, damping ratio and 
frequency of the first mode of vibration, respectively, and F and F p are the external and 
piezoelectric forces, respectively. The signal obtained from the accelerometer Vj, is 
proportional to the acceleration of the beam. 

v i “ aq (6.4.2) 

where a is the proportionality constant. A second order low pass filter is employed in the 
feedback loop which serves two purposes. First, high frequency noise, which could lead 
to unexpected results, is filtered out. Second, the low pass filter provides a 90° phase shift 
at its resonance, or cutoff frequency, so that the piezoelectric actuation forces are 180° out 
of phase with the displacement of the beam. Therefore, the control forces act as viscous 
dampers to reduce vibrational amplitudes in the beam. The governing equation for the 
second order low pass filter is as follows (Stout; 1976). 

v 0 + 2£, f to f v 0 + a>f v 0 = bv; (6.4.3) 

where vo, and cof are output voltage, damping ratio, and cutoff frequency of the filter and 
b is a constant of proportionality. The piezoelectric control forces are set to be proportional 
to the filter output voltage which is fed back into the control loop. 


(6.4.4) 



The block diagram for this system using the Laplace-transformed transfer functions for the 
above governing equations, assuming a=b=c=l, is presented in Fig. 6.1 1 The closed loop 
transfer function for this system is as follows. 


s + 2 £fCOfS + 00j- 


C(s) _ 

R(s) Ts + 2£ f co f s + )(s + 2^ b co b s + co b | + s 


(6.4.5) 


beam 



Fig. 6.1 1 Feedback control loop 

The control system is most effective when the cutoff frequency of the low pass filter is 
tuned to the natural frequency of the beam. 

“f"" 5 (6.4.5) 

A high filter damping ratio is also desirable to minimize the phase shift near the cutoff 
frequency. The damping ratio of the filter is selected to be that of a Butterworth second 
order low pass filter (£^-=0.7072) (Stout; 1976). The frequency response of the filter is 
presented in Fig 6.12 where the amplitude predicted by the transfer function is compared to 
experimental results with ojf = (25 Hz). The correlation of the filter response is excellent. 
The phase is also verified to be -90° at cof. 







Stability of the control system needs to be studied. This is especially true since the 
damping ratio corresponding to the structural damping in the beam for the first mode was 
found to be quite small (£f=0.003). A root locus of the control system is constructed to 
examine its stability characteristics with o f = co b = 157 rad/sec (25 Hz) (Fig. 6.13). This 
corresponds to parameters similar to beams 1-4. A similar analysis was also performed for 
beams 5-8. The two poles nearest the vertical axis correspond to the beam. They are very 
close the real axis due to the low damping in the beam. The other two poles correspond to 
the low pass filter. All values remain on the left hand side of the root locus plot indicating 
stability at any gain. This is a desirable feature of the control system selected to provide for 
increased damping. However, it is important to keep in mind that only a single bending 
mode has been used for the control design. It is possible that under certain circumstances, 
higher modes become excited at large gains which may threaten the stability of the control 
system. It is assumed for the current investigation that the influence of the higher modes is 
reduced by the low pass filter and do not present any stability problems. 




In practice, the natural frequency of the beam may not correspond exactly to the cutoff 
frequency of the filter. It is important to determine the effect of imprecise parameters on the 
control system. The root locus is presented in Figs 6.14 and 6.15 for beams frequencies 
which are 5% higher and 5% lower than the low pass filter cutoff frequency, respectively. 
If the open loop beam frequency is reduced by a small amount (Fig. 6.14), the poles 
corresponding to the beam are more strongly attracted to the single zero. The closed loop 
frequencies are calculated by finding the square root of the sum of the squares of the real 
and imaginary parts of the roots which are shown on the root locus. In this case, the root 
locus indicates that the closed loop frequency of the beam decreases as the gain is 
increased. Further analysis reveals that the closed loop damping increases to a value of 
approximately 0.3 until it starts to decrease at very high gains. When the beam natural 
frequency is slightly higher than that of the filter (Fig. 6.15), the root locus changes 
significantly. As before, the poles corresponding to the beam are nearest the vertical axis 
since the damping in the beam is small compared to the filter. However, these poles are no 
longer attracted to the single zero. Instead, they increase in magnitude until out of range of 
the figure. This indicates that the closed loop frequency of the beam increases, rather than 
decreases, as the gain is increased. As before, further analysis indicates that the closed 
loop damping increases to a value of approximately 0.3 until it starts to decrease at very 
high gains. 

Next, the frequency response of the control system is analyzed by constructing Bode 
plots of both the open and closed loop systems for the parameters considered above. (Figs. 
6.16 and 6.17). The open loop frequency response is presented in Fig. 6.16. A sharp 
peak exists in the amplitude response near the frequency of the beam. This is due to the 
low structural damping in the beam and was observed in the experimental open loop 
frequency response in the previous section. Only a single peak is observed here since only 
one mode is used in the control design. The phase shift near the natural frequency is very 
abrupt, again due to the low damping. The closed loop frequency response is presented in 



•• _ • . „ . : . ..144 

Fig. 6. 17. The amplitude peak near the natural frequency is not nearly as sharp as the one 

in the open loop Bode plot indicating increased damping. The magnitude of this peak is 
also reduced substantially. The phase shift near the natural frequency is much more 
gradual compared to the open loop response which also indicates increased damping due to 
the feedback control. 



Real Axis 


Fig. 6.13 Feedback control loop root locus 
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Fig. 6.16 Open loop frequency response. 



Fig. 6.17 Closed loop frequency response. 
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6.4.2 Stale space approach : Several modes must be included in the control analysis used 
for the experimental tests to more accurately determine the response of the coupled 
controls/structures system. A state space representation of the system is more convenient 
to work with than the classical approach and is developed in this section. However, the 
conventional notation of state space controls analysis used in Section 3.8 is not practical 
due to the presence of second derivative terms from the accelerometer output. Therefore, 
an alternative approach is used. The equations of motion of the structure are now given as 
follows. 

q + Z b q + A b q = <I> T (F + F p ) ( 6 . 4 . 6 ) 


where 




01 




(6.4.7) 


(6.4.8) 






(6.4.9) 


The matrix Ab contains the square of the natural frequencies, Zb contains the structural 
damping terms, <J> is the modal matrix containing the normalized eigenvectors and q is a 
vector of the modal paiticipation factors. The quantity F contains the applied forces and F p 
are the piezoelectric forces determined as follows. 


F p = GF p v 0 


(6.4.10) 
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where F p contains the piezoelectric forces per volt for one or more actuators with vq as the 

vector containing the applied voltages to the piezoelectric actuators and G contains the gains 
of each actuator. The signal obtained from the accelerometer Vj, is proportional to the 
acceleration of the beam as before. The signals from several accelerometers are related to 
the accelerations at different points as follows. 

Vj = SOq (6.4.11) 


where S contains the sensitivities of the accelerometers. The Endevco 2250A Micro 
miniature accelerometers used in this study have a sensitivity of 0.012 V/(m/s 2 ). The 
governing equation for a single second order low pass filter is also the same as before 
(Eqn. 6.4.3). Although only one filter is used in the current study, several filters can be 
used in a MIMO type control system with governing equations as follows. 

v 0 + Z f v 0 + A f v 0 - bVj (6.4.12) 


where the matrices Zf and Af are analogous to the matrices Zb and At>, respectively. The 
coupled controls/structures equations of motion are represented in state space as follows . 


I 

0 

0 

0 


q 


0 

I 

0 

0 

d 

q 


0 

0 

I 

0 

dt 

v o 


0 

-so 

0 

I 


.V 



i 

0 

0 • 

q 


0 

z b 

o t f p g 

0 

q 

+ 

o t f 

0 

0 

I 

v o 

0 

0 

-Af 

-Zf. 

_v 


0 


(6.4.13) 


The complex eigenvalues and eigenvectors of the above state space representation of the 
control system are determined from the following equation. 

(A - XB)x = 0 (6.4.14) 

where 

0 I 0 0) 

-A b -Z b O t F p G 0 

0 0 0 1 

0 0 -A f -Z f 


(6.4.15) 



10 0 0 
0 10 0 
0 0 10 
0 -SO 0 I 


J 49 

(6.4.16) 


Thus the state space representation used here reduces to an eigenvalue problem which is 
solved using standard techniques in a computationally efficient manner since no matrix 
inversions are necessary (Mignolet; 1995). Frequencies and damping ratios for the closed 
loop system are determined from the complex conjugate eigenvalue pairs as follows. 


<«„. = ^Re(y) 2 + Im(A. j )~ 


(6.4.17) 

(6.4.18) 


It is important to note that the above matrices, A and B, do not correspond to the A and B 
matrices often used in the state space approach for MIMO control systems. 


6.5 Closed Loop Structural Response 

In this section, the closed loop structural response is investigated. First, the 
experimental procedure is described and the design of the analog circuit to implement the 
control law are discussed. Results are then presented to correlate the predicted closed loop 
response using the higher order theory with experimental data. 


6.5.1 E xperimental procedure : The goal of the closed loop structural response investigation 
is to correlate the closed loop frequencies and damping ratios predicted by the HOT with 
experimental results. This was accomplished by recording the transient response of the 
beam due to an initial excitation. The hammer hit used in the open loop response to provide 
the initial excitation was not appropriate in the closed loop tests because this resulted in a 
large voltage spike at t = 0 which could have damaged the feedback control electronics and 
the data acquisition system. Instead, a random noise was used to excite the structure which 


was produced by applying a random voltage signal to piezoelectric transducer #2 for 
several seconds before the control loop was turned on (t < 0). At t = 0, the random noise 
was turned off and the feedback control loop was turned on. For t > 0, the output of the 
accelerometer was recorded for an appropriate length of time ( two seconds ) while the 
response of the system decayed. This procedure is illustrated in Fig. 6.18. The operation 
was repeated for all the test specimens for two different stacking sequences, four 
debonding lengths and seven different gains. Several tests were performed at different 
times to ensure consistency of the data. Several ARMAX models were generated for each 
test record to determine the frequencies and damping ratios. Engineering judgment was 
used to discard defective data. The remaining results for each test were then averaged. 
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Fig. 6.18 Test procedure. 


A5_.2 Analog circuit design : Implementation of the control law was a non-trivial procedure. 
The control loop needed to respond in real time to the motion of the structure while data 
was being taken at the same time Ambient noise was a factor and voltage and current limits 
of the actuators, data acquisition and control loop had to be considered. This lead to the 
implementation of the feedback control loop and test procedure via a nonlinear analog 
circuit which is shown in Fig. 6.19 (Durney et al; 1982). Not shown are the power 
supplies for the operational amplifiers and ground connections. Only the connections to the 
data acquisition system are shown for simplicity although the data acquisition software 
(Labview 3.0), hardware (Macintosh 7200/90) and I/O card (Labview PCI 1200) played an 
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integral role to the test procedure. The accelerometer also required a signal conditioner 

which is not shown. The entire circuit was designed, tested, and implemented in-house 
using commercially available electronic parts. The operation of the circuit is described 
next. 


B uflcr 



Fig. 6.19 Closed loop feedback circuit. 


At time t < 0, the data acquisition system outputs were as follows. 

DACOOUT = 1 volt, DACIOUT = 0 volt, (t < 0) 

These values positioned the relays in the circuit such that the random signal was connected 
to piezoelectric transducer #1 while the control loop remained disconnected. This provided 
the initial excitation to the structure and avoided a potentially dangerous voltage spike from 
the accelerometer. At time t = 0, these values swapped as follows. 

DACOOUT = 0 volt, DACIOUT = 1 volt ( 0 s; t <; tj ) 

This disconnected the random excitation and connected the feedback control loop. Data 
acquisition from the accelerometer output began immediately through ACHO and transpired 
for the required length of time (tj) to observe sufficient decay in the response of the 
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structure due to the feedback control. Once data acquisition was complete, these values 

returned to their original values, thus disconnecting the closed loop control, for safety. 

DACOOUT = 1 volt, DACIOUT = 0 volt ( 0 <; t <; t, ) 

The relays required 12 volts @ 1 1 ma to be triggered. However, the data acquisition 
outputs could only supply ± 5 volts @ 2 ma. Therefore, a separate power supply was used 
for the relay coils which are switched on and off using 2N2222A transistors The transistor 
bases were connected to both DACOOUT or DACIOUT through buffers to retain high 
impedance for protection. The relays were placed after the power amplifiers rather than 
before since the power amplifiers could potentially saturate from a floating input. 

The feedback control loop consisted of the accelerometer input, preamp, lowpass filter 
and power amplifier. The cutoff frequency of the low pass filter was adjusted for different 
test specimens by selecting appropriate values for the resistors and capacitors. It also 
provided the required phase shift at the resonant frequency of the beam. The gain of the 
preamp was selected by changing values of its feedback resistor. It provided a relatively 
high gain since the accelerometer signal level is so low (~10 ma). This high gain resulted 
in an annoying DC offset drift which was compensated for somewhat using another unity 
gain amplifier with an offset adjustment. 

The accelerometer output was also sent to the A/D converter in the data acquisition 
system (ACHO). It was amplified using a separate preamp from the feedback control loop. 
It was taken before the low pass filter since filtered data may have corrupt the system 
identification technique in post computation (Lyon; 1996). It was passed through a voltage 
clipper (± 5 volts) which utilized zener diodes to protect the data acquisition system. 

6.5.3 Transient response : Representative transient responses of beam 1 ( [0790°]3 S ,|3 = 0) 
are presented in Figs 6.20 (a-g) for increasing gains obtained from the experimental data. 
The quantity shown is the unfiltered output of the accelerometer which is in the millivolt 
range. Clearly, as the gain increases, the decay of the structure also increases. A smooth 
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exponential decay is observed for gains 0-2000 and the corresponding frequencies and 

damping ratios calculated by the ARMAX procedure were very consistent. At higher 
gains, the decay is not quite as smooth. Although the ARMAX frequencies were 
consistent, there was some variation in the damping ratios. This is due to the presence of 
nonlinear damping. The response of the analog circuit feedback control at higher gains 
may be partially responsible. The effect of higher modes may be enhanced at higher gains 
and also contribute to the appearance of the transient response at higher gains. It must be 
noted that the twist modes are neither controllable or observable in the current configuration 
and are assumed to have little effect on the transient response tests. 

The transient response decay envelopes for the first mode of vibration are presented in 
Figs. 6.21 (a-d) and 6.22 (a-d) for both [0790°]3 S and [457-45°]3 S beams. Debonding 
lengths of [3 = 0.0, 0.06, 0.12 and 0.18 for gains of 0, 200, 400, 1000, 2000, 4000 and 
8700 are shown. The response of each case is normalized to unity at t = 0 for comparison 
of the experimentally determined values using ARMAX (EXP) and the higher order theory 
(HOT). In general, agreement is very good. Increasing the gain significantly reduces the 
settling time for all cases. However, the presence of debonding tends to reduce the 
authority of the actuators and increases the settling time. 
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Fig. 6. 20 (a-d) Transient response of [0790°]3 S beam with increasing gain (p = 0). 
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Fig. 6.20 (e-g) Transient response of [0°/90°]3 S beam with increasing gain ([3 = 0). 






Mode 1 Mode 1 





6.5.4 Correlation of frequencies and damping ratios : In this section, the experimentally 
determined closed loop frequencies and damping ratios (EXP) are compared to the values 
of these quantities predicted by the higher order theory (HOT). Correlation is presented for 
both [0790 o ] 3 S and [45°/-45°]3 S beams with debonding lengths of [3 = 0.0, 0.06, 0.12 and 
0.18 for gains of 0, 200, 400, 1000, 2000, 4000 and 8700. The results are shown in 
Tables 6.6 and 6.7 and Figs 6.23 - 6.28. It must be noted that the mass of the 
accelerometer is included in the analysis using the HOT each case. 

Correlation between the closed loop frequencies for the EXP and HOT determined 
cases is very good in general as shown in Tables 6.6 and 6.7. Overall, increasing 
debonding length decreases the frequencies for both [0790°] 3s and [457-45°] 3s beams. 
This is because the debonded portion of the actuators does not contribute any stiffness to 
the overall structure while the mass is still present. Decreasing stiffness with constant mass 
results in reduced frequencies. This is also true since the debonded actuators have less 
control authority over the structure resulting in less influence of the control system on the 
structural response. The effects of increasing gain are observed most predominately in the 
nondebonded beams. The natural frequency of the [0790°] 3s beam ((3 = 0.0) increases 
11.1% from 25.2 to 28.0 Hz while the [457-45°]3 S (P - 0-0) beam increases 2.3% from 
17.3 to 17.7 Hz for the EXP case as the gains are increased from 0 to 8700. The HOT 
predicts similar increases in the natural frequency of 13.6% and 1.9%, respectively. The 
open loop frequencies of all of the [0790°]3 S test specimens are larger than the cutoff 
frequency of the low pass filter. Therefore, the closed loop frequencies of these beams are 
expected to increase with increasing gain as indicated by the root loci for the closed loop 
system presented earlier (Fig. 6.15). This trend is observed in both the EXP and HOT 
results. Although the open loop frequency of the nondebonded [457-45°]3 S beam is larger 
than the cutoff frequency of the low pass filter, the open loop beam frequency decreases 
significantly (~ 10%) for increasing debonding length. In fact, the open loop frequencies 
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of the debonded [457-45°] 3s beams (J3>0) are less than the cutoff frequency of the filter. 

Therefore, the closed loop frequencies are expected to decrease, rather than increase, 
according to the root loci presented earlier (Fig. 6.14). Although this effect is less 
noticeable for these beams since it is within the range of experimental accuracy, it is 
observed in the HOT results. 

The results for the closed loop damping ratios are presented in Figs. 6.23 and 6.24 for 
the [0790 °] 3s beams and Figs. 6.25 and 6.26 for the [457-45°] 3s beams with debonding 
lengths ranging from (3 = 0.0, 0.06, 0.12 and 0.18 and gains of 0, 200, 400, 1000, 2000, 
4000 and 8700 as before. Increases in closed loop damping ratios up to £ = 0.25 due to 
piezoelectric actuation are obtained. Agreement between the EXP and HOT results is 
excellent. The closed loop damping ratios include contributions from both the passive 
structural damping and the active control. The damping ratios corresponding to the passive 
structural damping using the HOT are set equal to the experimentally determined damping 
ratios at zero gain in each case These structural damping ratios are included in the 
calculation of the closed loop damping ratios. In general, debonding decreases the closed 
loop damping ratios, especially at high gain (Figs. 6.24 and 6.26). The presence of 
debonded actuators has the effect of increasing the open loop damping ratios (gain = 0) and 
low gain (gain = 200, 400) in some cases. This is due to friction between the interfaces of 
the debonded actuator and the substructure. The open loop damping ratios represent the 
passive structural damping present in the system. The structural damping dominates the 
closed loop response at low gain and explains why the damping may actually increase as 
the debonding length increases (Figs. 6.23 and 6.25). However, once the closed loop 
damping is above 0.01, the effect of the control system dominates the response and 
debonding consistently reduces the actuator authority in all cases. In the worst case, 
debonding reduces the closed loop damping ratio of the [0790°] 3s beam 68% from 0.023 to 
0.073 at a gain of 8700 as the debonding length is increased from (3 = 0.0 to {3 = 0.18. 
Debonding clearly has a detrimental effect on the closed loop damping ratios. This 
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detrimental effect on the control authority is especially clear on the settling time of the 

coupled controls/structures system as shown in Figs 6.27 and 6.28. The 2% settling time 
for the first mode is calculated as the time it takes for the control system to reduce the 


vibrational amplitude of the first mode to 2% of its initial value as follows. 
t _ ln(0.02) 


(6.5.1) 


where £b and cr>b are the closed loop frequencies corresponding to the first mode of the 
beam. It is presented in Figs. 6.27 and 6.28 for values of high gain at various debonding 
lengths for the [0°/90°]3 s and the [457-45°]3 S beams, respectively. Clearlv, increasing 
debonding significantly increases the settling time for both beams at relatively high gains. 
The settling time increases approximately 230% in the worse case for the [0790°] 3s beam 
as the debonding length increases from (3 = 0.0 to (3 = 0. 1 8 at a gain of 8700. 

It must be noted that, in general, the piezoelectric actuators have more control over a 
relatively soft structure compared to a relatively stiff one. The [457-45°] 3s beams are 
relatively softer than the [0790°]3 S beams in bending types of deformation. Therefore, the 
piezoelectric actuators should demonstrate greater control authority. However, the closed 
frequency is significantly less which reduces the output of the accelerometer. Therefore, in 
the current study, the feedback control demonstrates relatively less authority over the [457- 
45°]3 s beams. 



Table 6.6 First closed loop natural frequency of [0°/9Q°]3 S beams with debonding. 
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Table 6,7 First dosed loop natural frequency of [45%45°]3 S beams with debonding. 
Debonding Gain EXP HOT % error 

length ((3) 


CN 

VO 


’ — ' ' — 1 cN of co vo co 

^ b IT) IO IT) vd 


of 

vd 


CO 


of 

vd 


CO 

td 


of 

vd 


CO 


of of of in 

vd vd vd vd 


co co h 

d d d d 


o 

o o 

(N 


o 

o 

of 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

l> 

’ — 1 

CN 

of 

CO 


wo of of vo oo r-- cn \i- h vo oo \o « 

^ t >n cn cn cn cn oi cn cd 


vo vo 

WO WO* 


VO VO VO 

wd wd wd 


VO 

id 


VO 

id 


CN 

CN 

CN 

CN 

CN 

wd 

wd 

wd 

wd 

wd 


vd vo d vd \o vd \d 


VO VO vo vo VD VO 

wo wo wo wo wo wo 


vo 

wd 


o 

o o 

CN 


o 

o 

of 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 


’ — 1 

CN 

of 

CO 


o 


o 

o 

<N 


o 

o 

of 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

r- 

’ — 1 

CN 

of 

oo 


CN 'sj- CO CN of of 

d d d d d t d 


oo co oo oo oo oo C"~ 

vf d d ^ d ■d vt 


of 

of 

i/o 

of 

of 

of 

of 

wd 

wo 

wo 

wd 

wd 

wd 

wd 

i— 1 1 

i — i 

T— < 

< — i 

T—4 

r— 1 

i — i 


o 


o 

o 

CN 


o 

o 

of 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

r- 


CN 

of 

oo 


VO 

CN 

CO 

O 

i — i 


d 

o 

d 




Fig. 6.23 Closed loop damping ratios for Mode 1 of [0790°] 3s beams with debonding - low 

gain. 
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Fig. 6.24 Closed loop damping ratios for Mode 1 of [0790°]3 S beams with debonding - 

high gain. 
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Fig. 6.25 Closed loop damping ratios for Mode 1 of [45%45°]3 S beams with debonding 

low gain. 



Fig. 6.26 Closed loop damping ratios for Mode 1 of [45%45°]3 S beams with debonding 

high gain. 
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6,5.5 Laboratory setup : A visual record of the laboratory setup was obtained and is 

presented in this section. A test specimen is shown clamped in the vise ready to be tested 
in Fig. 6.29. Steel weights held the vise to a large concrete mass to minimize mechanical 
noise. The test specimen was mounted sideways to further reduce unwanted noise. The 
piezoelectric transducers are visible near the clamped end and the accelerometer can be seen 
near the tip of the beam. Some of the electronic devices used in the experiments are shown 
in Fig. 6.30. The two power amplifiers for the piezoelectric transducers are seen in the 
middle of the figure. On top of the power amplifiers are located the power supplies for the 
feedback control circuit. The feedback control circuit was constructed on two breadboards 
which are shown lying on the table. Connections between the power amplifiers, circuits, 
data acquisition system and test specimens were made via RG-58 cable with BNC 
connectors which comprise the “nest” of wires near the breadboards. 
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Fig. 6. 30 Feedback control electronic components 



7. Concluding Remarks 

A new theory has been developed to model composite laminates with surface bonded or 
embedded piezoelectric sensors and actuators. The theory is extended to incorporate the 
presence of pre-existing debonding in the composite laminate at the interface between the 
piezoelectric actuators and the underlying substructure. A refined higher order 
displacement field accurately captures the transverse shear deformation through the 
thickness of the smart composite laminate while satisfying the stress free boundary 
conditions on the free surfaces, including the debonded region. The higher order theory is 
implemented using the finite element method. The state space equations of motion utilizes 
the developed theory to account for piezoelectric sensing and actuation for closed loop 
control. Correlation of the higher order theory with analytical, numerical and experimental 
data is performed. Comparisons are also made with a commercial finite element code. An 
experimental investigation is conducted to further understand issues related to control of 
structures using piezoelectric actuation. Results are presented to demonstrate the analysis 
capability of the new theory. Additional studies are presented on the development of the 
constitutive equations for piezoelectric materials, justification of the higher order 
displacement field, and mesh generation which accounts for discrete piezoelectric layers. 
The following important observations are made. 

1) The developed higher order theory provides an accurate, computationally efficient 
analysis tool for the study of composite laminates with piezoelectric sensing and 
actuation. 


2) The theory is extended to incorporate debonded actuators. 
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3) The theory is implemented using the finite element method to allow incorporation of 

practical geometries, boundary conditions and discrete piezoelectric transducer 
locations. 

4) Correlation with analytical, numerical and experimental data is established for isotropic 
and orthotropic plates, piezoelectric actuation and debonding. Both static and dynamic 
results agree very well. 

5) The experimental investigation addressed practical issues, such as the design, 
construction, and testing of composite laminates and feedback control circuits. 
Consideration of these practical issues, coupled with the developed mathematical 
theory, provide an understanding and appreciation of the difficulty involved in 
implementation of smart structures as well as the potential benefit. 

6) The agreement between the higher order theory and the experiments is very good. 
Subtle trends in the natural frequencies with increasing gain and debonding length are 
experimentally correlated with the higher order theory. Closed loop damping ratios on 
the order of more than 20% of critical are obtained and agree with the developed theory. 
Debonding of the piezoelectric actuators results in a reduction of the control authority 
which is also correlated with the higher order theory. 

7) Significant deviations are observed between the higher order theory and the classical 
laminate theory in the strain and stress distributions through the thickness which 
indicates the importance of the higher order terms in the refined displacement field. 


. 

8) Debonding Of the piezoelectric layer, which is modeled using the higher order theory, is 
shown to cause peeling of the actuator away from the substructure. This has important 
implications for failure analysis. 


9) Local and global changes in the mode shapes also result due to debonding. 
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Appendix A 


Piezoelectric Constitutive Relations 



A.l General formulation 


Voigt (1928) formalized the equations for the free energy of a piezoelectric material 
which he called the “thermodynamic potential”. When the free energy is expressed in terms 
of strains, it is known as the first thermodynamic potential and is denoted by H(£,E). In 
tensor notation, it is formulated as follows (Tiersten; 1969). 

H(e,E) - ”Cij k [8ij£ ki — ejjj.Ej£jk — — TijEjEj (A.l) 

E c 

where Cy k i , ey k and rjjj are the coefficients of the elastic stiffness at constant electric field, 

piezoelectric stress and dielectric permittivity at constant strain, respectively. The quantities 
£;j and Ej are the components of strain and electric field, respectively. The first term on the 


right hand side of Eqn. A.l represent the elastic strain energy. The second term represents 
the electro-mechanical coupling while the third term is the electric energy. The differential 
coefficients of the first thermodynamic potential with respect to the components of elastic 
strain are the components of stress and the differential coefficients with respect to the 
components of electric field are the components of charge. These are the fundamental 
piezoelectric equations in terms of the strains and are given as follows 
3H(e,E) e i \ 

—^r = °ij = c ijki £ ki - e kij E k (converse effect) (A .2) 

ij 


dH(£,E) p / \ 

— — = Dj = e ikl £ kl + rj ik E k (direct effect) 


(A.3) 


The free energy can also be expressed in terms of stresses. It is then called the second 
thermodynamic potential, denoted by H(a,E) and is expressed as follows. 


H(o.E) = -I S 


ijklOijO 


kl +d ijk^i<7j k 



(A. 4) 


where s ki , dy k and r\f- are the coefficients of the elastic compliance at constant electric 
field, piezoelectric strain and dielectric permittivity at constant stress, respectively and Gy 
are the components of stress. The differential coefficients of the second thermodynamic 



potential with respect to the components of elastic stress are the conipoiients of strain and 
the differential coefficients with respect to the components of electric field are the 
components of charge These are the alternative form of the fundamental piezoelectric 
equations in terms of the stresses. 

3H(o,E) e / N 

3 a = e ij = s ijki a ki + d kij E k (converse effect) 


3H(g,E) 

5E; 


= Dj = d ikl a kl + T]? k E k (direct effect) 


(A .5) 
(A .6) 


The constitutive relations based on the first thermodynamic potential are used most often 

since they are a function of the components of strain rather than stress. However, the 
piezoelectric strain coefficients (d ikl ) are more commonly used than the stress coefficients 

( e iki ) • A relationship between the piezoelectric strain and stress coefficients allows the 

constitutive relations based on the first thermodynamic potential to be expressed using the 
piezoelectric strain coefficients. It is known that dielectric terms Tjfj and iq? are largely 

independent of stress and strain. Therefore, the following relationship using Eqns. A.3 
and A. 6 can be obtained. 

e ikl e kl “ d ikl a kl (A.7) 


The stress is related to the strain through Hooke’s law 
a ij ~ c ijkl 8 kl 


(A. 8) 


Substituting Hooke’s law into Eqn. A.7, the piezoelectric strain coefficients are now 
expressed as follows. 

e ikl = d iki c klmn (A. 9) 


Substituting the above equation into Eqns. A.2 and A.3 yields the following constitutive 
relations expressed in terms of the components of strain and piezoelectric strain 
coefficients. 
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CT ij ~ c i.'kl c ki dkijConmEv (converse effect) (A. 10) 

D i = d ikl c kimn e mn +1 4 E k (direct effect) (A. 11) 

Rearranging and using matrix notation produces the constitutive equations for piezoelectric 
materials. 

a = Q^e - d T E) (converse effect) (A. 12) 

D = TE + dQe (direct effect) (A. 13) 

where Q is the elastic stiffness matrix, d is the piezoelectric strain coefficient matrix, Y is 
the dielectric constant matrix, a and 8 are the stress and strain vectors, respectively and D 
and E are the electric charge and the electric field vectors, respectively. 

A.2 Theoretical background 

The focus of this dissertation is to utilize the macroscopic properties of piezoelectric 
materials and integrate them as elements of smart composite structures as both sensors and 
actuators. To achieve this goal, it is necessary to understand the electromechanical 
constitutive relationships which govern these materials . As mentioned previously, Jacques 
and Pierre Curie discovered the direct piezoelectric effect in 1880 when they observed that 
an electric charge developed on certain materials in response to a mechanical stress. Soon 
thereafter, Gabriel Lippman predicted the converse piezoelectric effect which occurs when a 
mechanical strain develops in response to an applied electric charge. His rational was as 
follows. Consider the following thermodynamic potential function formulated in terms of a 
single scalar stress and electric field (Cady; 1964). 

C - — c a 2 + -r|'E 2 + d'oE 
2 2 


(A. 14) 


where c is the elastic compliance, r\' is the permittivity, a is the stress, E i s the el eel ric fi eld 
and d' is the piezoelectric strain coefficient. The change in £ due to a small variation in E 
and c is as follows. 


d^ = DdE + eda 


where D is the electric charge and £ is the strain which are found as follows. 


D = 


3E 

K 

3a i 


(A. 15) 


(A. 16) 
(A. 17) 


Since the order in which differentiation occurs, it follows that 

3D _ 3e 
3a E 3E a 


(A. 18) 


Since it is known that the relationship between an applied stress and an electric charge is 
linear over a wide range of pressure for most crystals, the following is true. 



(A. 19) 


where 5 is a constant. This was verified by the Currie brothers to be true for piezoelectric 
crystals as well. Lippman showed that not only must converse relationship exist between 
the strain and an applied electric field, it must be equal to the same constant previously 
found for the direct relationship. 




(A .20) 


His theoretical prediction was soon found to be true through further experiments by the 
Curries. It must be noted that Eqn. A. 15 is analogous to the enthalpy of a reversible 
system where the quantities D, E, £ and a are analogous to absolute temperature, entropy, 
volume and pressure in ordinary thermodynamics. The electric charge, D, is also 
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analogous to the strain, e, and is Often termed the electric displacement. Similarly, the 
electric field, E, is analogous to the stress, a, and is often termed the electric force. The 
above analysis can easily be expanded to incorporate a multi dimensional state of stress and 
electric field. 

A.3 Physical interpretation of piezoelectric constants 

Several coefficients are defined to characterize the properties of a piezoelectric material. 
Three axes are used to identify directions in a piezoelectric material which are necessary to 
define these coefficients. These axes, termed 1 , 2 and 3, are analogous the X, Y and Z of 
the classical three dimensional orthogonal set of axes and are indicated in Fig. A.l (a). The 
polar, or 3, axis is taken parallel to the direction of polarization within the ceramic. This 
direction is established during manufacturing by a high DC voltage that is applied between 
a pair of electroded faces to activate the material. The polarization vector “P” is represented 
by an arrow pointing from the positive to the negative poling electrode. Piezoelectric 
coefficients with double subscripts link electrical and mechanical quantities. The first 
subscript gives the direction of the electrical field associated with the voltage applied or the 
charge produced. The second subscript gives the direction of the mechanical stress or 
strain. For the current research, the most important coefficients are the piezoelectric strain 
coefficients (dy) which relate the mechanical strain developed in a piezoelectric material to 

the applied electric field. 

, strain developed 

applied electric field J 

Conversely, the dy coefficients also relate the charge collected due to an applied mechanical 
stress. These are contained in the piezoelectric strain coefficient matrix, d, which has the 
following form for piezoelectric material having orthorhombic mm2 symmetry such as 
PVDF (polyvinylidene fluoride) and PZT (lead zirconate titanate). 



0 0 d 31 

0 0 d 32 

0 0 d 33 

0 d 24 0 

d i5 0 0 

0 0 d 36 _ 
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(A .22) 


The dimensions are expressed as length per length, per volts per length (length per volt). 
For instance, the d 3 i coefficient relates the coupling between an applied electric field 
(voltage) in the 3 direction to the strain in the 1 direction as shown in Fig. A.l (b). 



( a ) (b) 

Fig. A1 . Piezoelectric material: (a) coordinate system (b) effect of d 3 j coefficient. 


The piezoelectric stress coefficients, gy, relate the 

mechanical stress as follows. 
o _ open circuit electric field 
° applied mechanical stress 

The dimensions for gy are expressed as volts per length, per force per length squared. 
These coefficients are used less often then the dy coefficients since normally the electric 
field is taken as the independent variable. However, their is a direct correspondence of 
nonzero values of dy to nonzero values of gy. The relative dielectric constant K is the 
dimensionless ratio of the permittivity of the materials, to the permittivity of free space, 


open circuit electric field to the applied 


(A .23) 




The piezoelectric strain and stress coefficients are related through the dielectric properties as 
follows. 



Appendix B 


Displacement Field Justification 
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-The formulation of the higher order displacement field is based on both physical 

relevance and considerations for the finite element implementation in order to solve useful 
problems. A great deal of intuition is gained by solving a simple beam problem, using the 
Ritz method, which is used to illustrate concepts which can be applied to the advanced 
formulation of the smart composite laminate. The beam has length L, height h, distributed 
load q and elastic and shear moduli E and G, respectively. One end is fixed as indicated in 
Fig. B.l. 


q 


p 1 

a 


z 

z 

z 

z 

Z 

Z 


/ 

A 





E,G 



kd 




Fig. B 1 Cantilever beam with distributed load. 


The general higher order displacement field for the beam is given as a cubic expansion 
of the inplane displacements (u) and a single function for the transverse displacement (w) 
as follows. 

U = 4>io( x ) + 1 (x) + Z 2 C() j2 (x) + 2 3 c|)i 3 (x) 

i / x (B.l) 

W = <ho(x) 


The in plane term (j> 10 (x,y) is ignored for this simple investigation. The relevant normal 
and shear strain are given as follows. 


3u 



du dw 
dz dx 


(B.3) 


where 8 is the normal strain and T is the transverse strain. The total potential energy is 
given as follows. 




Lb h/2 

— n = -JJ J [lie 2 + Gx 2 jdzdydx - qw 

^00 -h/2 
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(B.4) 


The principle of minimum total potential energy states that the equilibrium position is found 
which minimizes the potential functional Id. This is achieved by applying standard 
procedures in variational calculus which are used to find a stationary point of II and derive 
the Euler-Lagrange equations which constitute the equilibrium equations and boundary 
conditions. The results are well known and not presented here for the sake of brevity. 
Instead, approximate solutions are found using the Ritz approach to determine values for 
the coefficients of interpolating polynomials. These approximate solutions represent 
properties of the finite element solution which is used for implementation of the higher 
order theory (Chapter 3). This is achieved by finding a stationary point to II with respect 
to the unknown coefficients of the approximation polynomials (aj) as follows. 


an 

da.; 


-0 


(B.5) 


Using this approach , the coefficients can be determined by solving the set of resulting 
equations from Eqn. B.5. 


Classical theory 

First, consider the case of classical theory which can be obtained by setting 
§\i( x ) = 0i3(x) = 0. The quantity (j) 12 (x,y) is determined by using the Kirchcoff 

condition that states planes normal to the neutral axis remain plane after bending for 
sufficiently thin beams and the shear stress vanishes as follows. 

T=4 ’> 1+ ^ i=0 (B - 6) 


Setting (J) 31 = w 0 yields the following displacement field for classical beam theory. 




W = Wg 
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(B.7) 


Selecting a two term Ritz solution which satisfies the geometric boundary conditions at the 

2 3 

fixed end, w 0 = b 2 x. + b 3 x , results in the following values for the unknown 


interpolation polynomial. 


w n - 


5qL^ 

2 Eh 2 


_qL 3 
Eh 2 


(B.8) 


A three term Ritz solution, w 0 = b 2 x“ + b 3 x 3 + b 4 x 4 , yields the exact solution. That is, 


w 0 = 3^x 2 -2-5t-x 3 +--^x 4 


Eh' 


Eh' 


2 Eh' 


(B.9) 


Upon closer examination, the two term solution is a least squares fit to the exact solution. 
It gives the exact solution at the points where the second order Legendre polynomial 
vanishes. It is important to note that the shear stiffness, G, does not enter here since the 
shear strain is assumed to be zero. The theory is valid for thin beams where the shear 
strain is not important. However, significant errors are obtained for thicker beams where 
the shear strain becomes important. 


First order theory 

The first order theory, which is identical to the Timeshenko beam theory, accounts for 
additional deformation of the cross section. It recognizes that the action of the shear force 
causes a shear strain which in turn causes a warping of planes normal to the neutral axis. 
The theory in effect averages the effect of shear strain over the entire cross section. This 
theory utilizes two functions, (j) n and (}) 31 . Setting <j>]j =u, and <j) 30 = w 0 yields the 

following displacement field for first order theory. 

U = ZUj 
W = Wq 


(B. 10) 
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Selecting a two term Ritz solution which satisfies the geometric boundary conditions at the 

fixed end, u, - ajX and w 0 = bjX , now results in the following values for the unknown 
interpolation polynomial. 


Un = -3- 


qL^ 


w 


0 GL 2 +Eh 2 

1 qL(4GL 2 + Eh 2 ) 


0 


2 G(GL 2 + Eh 2 ) 


(B.ll) 


(B.12) 


A thin beam in which the shear strain is not important can be represented by an infinite 
shear stiffness (G -» «>). It is clear that as G -> <*>, both a! and bj -» 0. This serious 
deficiency in the solution and a clear case of a shear “locking” since a large load q produces 
practically no displacements for thin beams. The principle of minimum potential energy 
can be viewed as a minimization process of the strain energy potential. The shear term in 
the energy functional, Gx , represents a constraint, x, with a penalty term, G. As the limit 
of this penalty term approaches infinity, a penalty condition is introduced on the shear 
strain which requires that the shear strain x = a 1 x + b i must vanish in the Ritz 
approximation which is the case for thin beams. This penalty term requires that both aj and 
bi — > 0 which also forces the bending strain to zero and does not accurately represent the 
physical situation. A meaningful solution can not be obtained in this case which results 
from the inconsistent interpolation of the shear strain. A three term Ritz solution with 
U| = ajx and w 0 = bjX + b 2 x , gives the following approximate solution. 


0 qkr 

U n = —2 ~ 7T X 

Eh 2 


qL l q(-2GL 2 + Eh 2 ) 


w 0 = — X 

u G 2 


EGh 2 


(B. 1 3) 


(B. 14) 


The shear strain is now given as 
x = bj + (&] +2b 9 )x 


(B.15) 
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N o w ,• as- G -^> oo the shear strain penalty term governed by G also — > 0 since bi — > 0 arid 
(aj + 2b 2 ) — » Oas well. Therefore, the shear strain is now effectively approximated, even 

for a thin beam with a rigid shear stiffness. This is because the shear strain does not 
introduce a spurious constraint so that the bending strain, although constant, is not reduced 
to zero and the physical aspect of the problem is accurately modeled. The reason this 
interpolation is effective is because the constraint (aj + 2b 2 ) is consistently balanced. That 

0Wn 

is, both u 0 and are of the same order. A finite element model based on this 

dx 

interpolation scheme is completely free from locking. However, a finite element model 
which uses quadratic interpolation is quite cumbersome to implement since it requires that 
extra nodes be defined which are not required using either linear or cubic interpolations. 
Another drawback is the representation of the shear strain as a linear function. Shear 
correction factors are required to obtain an accurate solution since the shear is a linear 
approximation of a quadratic function. The shear strain distribution also violates the stress 
free boundary conditions on the free surfaces which require the shear strain to vanish on 
the top and bottom surfaces of the beam. 


Higher order theon 


The conditions that the shear strains vanish are used by the higher order theory to 
satisfy this requirement. Consider the following displacement field, 
u = z<t>, !(x) + z 2 <t> 12 (x) + z 3 <j> l3 (x) 

W = (j> 3 0 (x) 


In this general form, this higher order displacement field does not satisfy the requirement 
that the shear stresses vanish on the free surfaces. Imposing this requirement is equivalent 
to setting the shear strain equal to zero at ± h/2. 

£ lh/2 = £ L/2 = 0 (B.17) 


This results in the following equations 



12) :: ,Q,.= ii i(x) + h<j), 2 (x) + h 2 (J> ! 3 (x) + 

4 dx 
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(B.18) 


T(-h/2) = 0 = (|) ll (x)-h<ti 12 (x) + ^h 2 (|) l 3(x) + ^flW ( B . 19 ) 

4 dx 

Clearly, <t> 12 (x) must be zero to satisfy both equalities. In addition, the higher order 
function tj>] 3 (x) can be solved for in terms of the lower order functions. This leads to the 


following refined displacement field. 


4z 3 ( 


u = Z<l>ll(x)- 2 
3h 

w = <t>3o( x ) 


<j) n (x) + 


^30 (x) 

dx 


(B.20) 


The shear stress is then a quadratic function of z given by the following 


T = 


1-4- 


.2 A/ 




<t>l]W + 


dx 


(B.21 ) 


which vanishes at ±h/2 as required. The refined theory, which contains the same number 
of unknown functions as first order theory, is finally obtained by setting (j ) 15 — iq and 
$30 ~ w o as follows. 


u = ZU](x) - 
W - Wq(x) 



dwo(x)^ 

dx J 


(B.22) 


However, the appearance of the derivative term in u requires a higher order interpolation of 
w 0 . Choosing a three term Ritz solution, Uj =ajx and w 0 = b 2 x 2 + b 3 x 3 , yields the 

following approximate solution. 

_ 10qL 2 (-7056G 2 L 4 - 2261GEL 2 h 2 + 100E 2 h 4 ) 

Eh 2 (35280G 2 L 4 + 7336GEL 2 h 2 + 25E 2 h 4 ) 


(B.23) 



5q!. 5 (l 4 ) iiG 2 ']; 1 -f 49.S0 iGI-:i. ;! h 2 + 11 65E 2 h 4 ) J 
2Eh 2 (35280G 2 L 4 + 7336GEL 2 h 2 +25E 2 h 4 ) 

(l68GL 2 + Eh 2 ) 

525qL 7 — L — jy x 2 

(35280G 2 L 4 + 7336GEL 2 h 2 + 25E 2 h 4 ) 



(B.24) 


Taking the limit G — > oo to represent a thin beam which has rigid shear stiffness of the Ritz 


solution yields the following displacement field. 

u =-2-Sk-x 
Eh 2 


w = 


Eh" 


(B.25) 

(B.26) 


Despite the additional complexity of the refined displacement field, this approach reduces to 
the three term Ritz solution for first order theory. It is also no better than a two term 
solution for the refined displacement field using only the quadratic term for w since b 2 
— » 0 for thin beams. In the finite element implementation of the refined displacement 
field, the quadratic interpolation of w is cumbersome to implement. A cubic interpolation 
that introduces additional degrees of freedom is not justified either since no improvement is 
observed in the solution over the simpler first order theory for thin plates. A better 
approach is needed. 


Modified higher order theory 


Considering the physical nature of the 

( d w 

result from the rotation of the beam 


3x 


problem, it is expected that the bending strain 
\ plus a distortion effect due to shear Uj , plus 


some higher order terms which can be solved in terms of lower order terms by using the 

stress free boundary conditions. This situation can be modeled mathematically by setting 
3wa 

0U =Uj — - — and (j) 30 = w 0 . Substituting into the refined displacement field (Eqn. 
ox 

B.22) yields the following modified refined displacement field 



u, - 


3wo(x)V -1z 


3x 


3h 


:2 


u l(x) 


W = W 0 (x) 
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(B.27) 


Using the same three term approximation to the unknown functions, u ] =ajX and 
w o - ^ 2 X + t*3 x ’ yields the following approximate solution for the modified refined 


approach. 

210qL 2 

u= -■■■ x 

(280GL 2 +Eh 2 ) 


(B.28) 


qL 2 (l400GL 2 + 173Eh 2 ) 2 qL 3 
2Eh 2 (280GL 2 +Eh 2 ) X Eh 2 X 


(B.29) 


Again, taking the limit G °° to represent a thin beam and reduce the shear constraint to 
zero yields the following 


u - 0 


Wr 


5 q L£ 
2 Eh 2 


q L ..3 
Eh 2 


(B.30) 
(B .3 1 ) 


These equations are exactly the same as the three term Ritz solution using the classical 
theory where the shear stress is assumed to be zero which is true for a beam with infinite 
shear stiffness. Therefore, the modified displacement field yields the correct trend for thin 
beams by reducing the shear constrain to zero. The higher order interpolation polynomials, 
which require additional degrees of freedom in the finite element interpolation, are now 
justified since the bending rotation is now a linear function rather than a constant as in the 
first order theory. The shear strain is given as follows 


T = Ui - 


dw 0 i ^w 0 


(B.32) 


which, when calculated in terms of the Ritz coefficients, t = aj, vanishes for thin plates as 
expected. For thicker plates, the shear strain is accurately accounted for as a quadratic 
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function which vanishes on the frce^surfaces;' Consistent interpolation for the shear strain" 

is automatically satisfied while avoiding quadratic interpolation schemes which are difficult 

to implement. Understanding the physical nature of the problem and viewing the principle 

of minimum potential energy as a minimization process which includes the shear as a 

constraint allows a modified refined displacement field to be formulated. This displacement 

field accurately accounts for the shear strain and produces no spurious constraint which 

frees the solution from any form of locking. 


Appendix C 


Mesh Generation 
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Each element in the finite element model forthe current study is formulated assuming 

that the piezoelectric layers are either present throughout the entire element, or are 
completely absent. While it is possible to formulate elements which only partially contain 
piezoelectric layers, it is quite laborious and computationally inefficient. Since the finite 
elements either contain, or do not contain, piezoelectric layers, the boundaries of 
piezoelectric layers must coincide with the element nodes. Therefore, a change in the shape 
of a piezoelectric layer must be reflected by a change in the finite element mesh. This task 
is non trivial and impractical to accomplish by hand for realistic mesh sizes or at each 
iteration of a closed loop optimization procedure. Therefore, an efficient strategy must be 
devised to generate a finite element mesh to accommodate arbitrary piezoelectric layer 
geometries. 

The first step in accomplishing the mesh generation is to choose a proper mesh for the 
composite laminate with no piezoelectric layers. Ideally, this mesh should have equally 
sized elements for simplicity. The aspect ratio of the k-th element is defined as follows, 
lx 

a k=7r~ (C.l) 

V 

where and ly k are the length of the k-th element in the x and y direction, respectively. 
These geometric parameters are calculated from the nodal coordinates of each element as 


follows (Fig. C.l). 


■x t = x j+1 - X, 

(C.2) 

ki = y j+i - yj 

(C.3) 


The element aspect ratio should be as close to one as possible since the accuracy of the 
finite element solution degrades as the aspect ratio differs from one. However, this equally 
sized mesh does not necessarily coincide with the boundaries of an arbitrary piezoelectric 
layer. Consider the introduction of a piezoelectric layer on a square laminate as shown in 
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•' Fig. C.l (a) The original mesh of square elements is shown in Fig. C. 1 (b) along with the 

modified mesh to incorporate the piezoelectric layer. In this case, it is clear that node points 

of the original mesh corresponding to X 3 and y 2 must be changed to accommodate the 

piezoelectric layer. In meshes with a practical number of elements, the process of 

determining nodal coordinates to incorporate the piezoelectric layer while retaining an 

aspect ratio as near to one as possible is a non trivial task. 



Fig. Cl Finite element mesh incorporating piezoelectric layers. 


The process of determining the best possible mesh for a given actuator size and location 
is an optimization problem within itself. The objective is to retain the original element 
aspect ratios as closely as possible of the original mesh while incorporating the 
piezoelectric layer (an aspect ratio of one is best). A quadratic penalty function is defined 
which the sum of the squares of the difference between the element dimensions of the old 


and new mesh as follows 
NELEM , 

1 


=7 X 


I e 

U J *k 


2 +(] e ' -l e ' 2 
l Yk Yk 


k = l 


(C.4) 


The quantities 1* and 1® are the element aspect ratios for the new mesh which are 
calculated as follows. 


i e — y' — y' 
1 x k A i+1 A i 


(C.5) 



\ 
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(C.6) 


The quantities x( +1 , xj, yj +1 and yj are the nodal coordinates for the new mesh. 

Minimizing this penalty function corresponds to retaining the original aspect ratios as 
closely as possible. Of course, constraints need to be imposed to ensure that the original 
geometry of the laminate is retained and that the proper nodal values correspond to the 
edges of the piezoelectric layer. These constraints will be discussed shortly. The minimum 
value of F is found by setting the derivatives with respect to the new nodal coordinates to 
zero. First, F is rewritten in terms of the nodal coordinates for the original and new mesh. 


1 m-in . 


F = 2 XX (W+i " _ fa+i ~ x il) + ({yj+i {yj+i - y } }] 

i-l j-1 


(C.7) 


where M and N are the number of nodes in the x and y directions, respectively. Next, 
derivatives of F are calculated with respect to the new nodal coordinates and are set to zero 


as follows. 




3F 

— = -x' +1 +2x'-x;_ 1 =0 

i = 2,3,- 

•*,M - 1 

(C.8) 

3y<- yj+i +2 u y h = 0 

j = 2,3,- 

•*,N - 1 

(C.9) 


where it is assumed that all of the nodes in the original mesh are evenly spaced in the x and 
y directions. It must be noted that the coordinates of the edges of the laminate are 
considered fixed (i=l,M and j=l,N) to retain the original geometry of the laminate. 
Therefore, derivatives of the penalty function with respect to the nodal coordinates which 
correspond to these nodes are zero. Four additional constraint equations are required to 
ensure that this is true. 

xj = 0, x' M =L (C.10) 


yl=0, y' N =W 


(C.ll) 
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where L and W are the length and width of the laminate, respectively. Constraints are also 
imposed on the nodes nearest the edges of the piezoelectric layer to ensure that these nodes 
are aligned with the border of the piezoelectric layer. This is accomplished by setting to 
zero all coefficients of the relevant equation except for a single term corresponding to the 
constrained node. The term on the right hand side is then set to the value of the nodal 
coordinate corresponding to the constrained value of the nodal coordinate. The above 
linear equations are uncoupled in the sense that the nodal coordinates in the x direction do 
not depend on the nodal coordinates in the y direction and vice versa. The equations which 

corresponds to the nodal coordinates in the x direction are written in matrix form as follows 

" 1 

-1 2-1 O' s 

1 

-1 2 -1 

1 

0’s -12 -1 

1 



x ! 


0 


x 2 


0 


x 3 


a l 

< 

x 4 

> = 4 

0 

> 


X M~2 


a 2 


X M-1 


0 


X M . 


L 


(C.12) 


where the constraint equations are flagged with an asterisk (*). Although the resulting set 
of equations are not symmetric, they can easily be transformed into a symmetric set of 
equations by performing a few elementary row operations. Solving this set of equations 
yields optimal values of the nodal coordinates in the x direction to retain element aspect 
ratios as near to the original values as possible while incorporating new coordinates due to 
the introduction a piezoelectric layer. The equations to determine the nodal coordinates in 
the y direction are similar. 



